home *** CD-ROM | disk | FTP | other *** search
/ CU Amiga Super CD-ROM 6 / CU Amiga Magazine's Super CD-ROM 06 (1996)(EMAP Images)(GB)(Track 1 of 4)[!][issue 1997-01].iso / cucd / prog / gnu-c / src / gcc-2.7.0-amiga / config / m68k / lb1sf68.asm < prev    next >
Assembly Source File  |  1995-06-15  |  81KB  |  2,977 lines
  1. /* libgcc1 routines for 68000 w/o floating-point hardware. */
  2. /* Copyright (C) 1994 Free Software Foundation, Inc.
  3.  
  4. This file is free software; you can redistribute it and/or modify it
  5. under the terms of the GNU General Public License as published by the
  6. Free Software Foundation; either version 2, or (at your option) any
  7. later version.
  8.  
  9. In addition to the permissions in the GNU General Public License, the
  10. Free Software Foundation gives you unlimited permission to link the
  11. compiled version of this file with other programs, and to distribute
  12. those programs without any restriction coming from the use of this
  13. file.  (The General Public License restrictions do apply in other
  14. respects; for example, they cover modification of the file, and
  15. distribution when not linked into another program.)
  16.  
  17. This file is distributed in the hope that it will be useful, but
  18. WITHOUT ANY WARRANTY; without even the implied warranty of
  19. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20. General Public License for more details.
  21.  
  22. You should have received a copy of the GNU General Public License
  23. along with this program; see the file COPYING.  If not, write to
  24. the Free Software Foundation, 59 Temple Place - Suite 330,
  25. Boston, MA 02111-1307, USA.  */
  26.  
  27. /* As a special exception, if you link this library with files
  28.    compiled with GCC to produce an executable, this does not cause
  29.    the resulting executable to be covered by the GNU General Public License.
  30.    This exception does not however invalidate any other reasons why
  31.    the executable file might be covered by the GNU General Public License.  */
  32.  
  33. /* Use this one for any 680x0; assumes no floating point hardware.
  34.    The trailing " '" appearing on some lines is for ANSI preprocessors.  Yuk.
  35.    Some of this code comes from MINIX, via the folks at ericsson.
  36.    D. V. Henkel-Wallace (gumby@cygnus.com) Fete Bastille, 1992
  37. */
  38.  
  39. /* These are predefined by new versions of GNU cpp.  */
  40.  
  41. #ifndef __USER_LABEL_PREFIX__
  42. #define __USER_LABEL_PREFIX__ _
  43. #endif
  44.  
  45. #ifndef __REGISTER_PREFIX__
  46. #define __REGISTER_PREFIX__
  47. #endif
  48.  
  49. #ifndef __IMMEDIATE_PREFIX__
  50. #define __IMMEDIATE_PREFIX__ #
  51. #endif
  52.  
  53. /* ANSI concatenation macros.  */
  54.  
  55. #define CONCAT1(a, b) CONCAT2(a, b)
  56. #define CONCAT2(a, b) a ## b
  57.  
  58. /* Use the right prefix for global labels.  */
  59.  
  60. #define SYM(x) CONCAT1 (__USER_LABEL_PREFIX__, x)
  61.  
  62. /* Use the right prefix for registers.  */
  63.  
  64. #define REG(x) CONCAT1 (__REGISTER_PREFIX__, x)
  65.  
  66. /* Use the right prefix for immediate values.  */
  67.  
  68. #define IMM(x) CONCAT1 (__IMMEDIATE_PREFIX__, x)
  69.  
  70. #define d0 REG (d0)
  71. #define d1 REG (d1)
  72. #define d2 REG (d2)
  73. #define d3 REG (d3)
  74. #define d4 REG (d4)
  75. #define d5 REG (d5)
  76. #define d6 REG (d6)
  77. #define d7 REG (d7)
  78. #define a0 REG (a0)
  79. #define a1 REG (a1)
  80. #define a2 REG (a2)
  81. #define a3 REG (a3)
  82. #define a4 REG (a4)
  83. #define a5 REG (a5)
  84. #define a6 REG (a6)
  85. #define fp REG (fp)
  86. #define sp REG (sp)
  87.  
  88. #ifdef L_floatex
  89.  
  90. | This is an attempt at a decent floating point (single, double and 
  91. | extended double) code for the GNU C compiler. It should be easy to
  92. | adapt to other compilers (but beware of the local labels!).
  93.  
  94. | Starting date: 21 October, 1990
  95.  
  96. | It is convenient to introduce the notation (s,e,f) for a floating point
  97. | number, where s=sign, e=exponent, f=fraction. We will call a floating
  98. | point number fpn to abbreviate, independently of the precision.
  99. | Let MAX_EXP be in each case the maximum exponent (255 for floats, 1023 
  100. | for doubles and 16383 for long doubles). We then have the following 
  101. | different cases:
  102. |  1. Normalized fpns have 0 < e < MAX_EXP. They correspond to 
  103. |     (-1)^s x 1.f x 2^(e-bias-1).
  104. |  2. Denormalized fpns have e=0. They correspond to numbers of the form
  105. |     (-1)^s x 0.f x 2^(-bias).
  106. |  3. +/-INFINITY have e=MAX_EXP, f=0.
  107. |  4. Quiet NaN (Not a Number) have all bits set.
  108. |  5. Signaling NaN (Not a Number) have s=0, e=MAX_EXP, f=1.
  109.  
  110. |=============================================================================
  111. |                                  exceptions
  112. |=============================================================================
  113.  
  114. | This is the floating point condition code register (_fpCCR):
  115. |
  116. | struct {
  117. |   short _exception_bits;    
  118. |   short _trap_enable_bits;    
  119. |   short _sticky_bits;
  120. |   short _rounding_mode;
  121. |   short _format;
  122. |   short _last_operation;
  123. |   union {
  124. |     float sf;
  125. |     double df;
  126. |   } _operand1;
  127. |   union {
  128. |     float sf;
  129. |     double df;
  130. |   } _operand2;
  131. | } _fpCCR;
  132.  
  133.     .data
  134.     .even
  135.  
  136.     .globl    SYM (_fpCCR)
  137.     
  138. SYM (_fpCCR):
  139. __exception_bits:
  140.     .word    0
  141. __trap_enable_bits:
  142.     .word    0
  143. __sticky_bits:
  144.     .word    0
  145. __rounding_mode:
  146.     .word    ROUND_TO_NEAREST
  147. __format:
  148.     .word    NIL
  149. __last_operation:
  150.     .word    NOOP
  151. __operand1:
  152.     .long    0
  153.     .long    0
  154. __operand2:
  155.     .long     0
  156.     .long    0
  157.  
  158. | Offsets:
  159. EBITS  = __exception_bits - SYM (_fpCCR)
  160. TRAPE  = __trap_enable_bits - SYM (_fpCCR)
  161. STICK  = __sticky_bits - SYM (_fpCCR)
  162. ROUND  = __rounding_mode - SYM (_fpCCR)
  163. FORMT  = __format - SYM (_fpCCR)
  164. LASTO  = __last_operation - SYM (_fpCCR)
  165. OPER1  = __operand1 - SYM (_fpCCR)
  166. OPER2  = __operand2 - SYM (_fpCCR)
  167.  
  168. | The following exception types are supported:
  169. INEXACT_RESULT         = 0x0001
  170. UNDERFLOW         = 0x0002
  171. OVERFLOW         = 0x0004
  172. DIVIDE_BY_ZERO         = 0x0008
  173. INVALID_OPERATION     = 0x0010
  174.  
  175. | The allowed rounding modes are:
  176. UNKNOWN           = -1
  177. ROUND_TO_NEAREST  = 0 | round result to nearest representable value
  178. ROUND_TO_ZERO     = 1 | round result towards zero
  179. ROUND_TO_PLUS     = 2 | round result towards plus infinity
  180. ROUND_TO_MINUS    = 3 | round result towards minus infinity
  181.  
  182. | The allowed values of format are:
  183. NIL          = 0
  184. SINGLE_FLOAT = 1
  185. DOUBLE_FLOAT = 2
  186. LONG_FLOAT   = 3
  187.  
  188. | The allowed values for the last operation are:
  189. NOOP         = 0
  190. ADD          = 1
  191. MULTIPLY     = 2
  192. DIVIDE       = 3
  193. NEGATE       = 4
  194. COMPARE      = 5
  195. EXTENDSFDF   = 6
  196. TRUNCDFSF    = 7
  197.  
  198. |=============================================================================
  199. |                           __clear_sticky_bits
  200. |=============================================================================
  201.  
  202. | The sticky bits are normally not cleared (thus the name), whereas the 
  203. | exception type and exception value reflect the last computation. 
  204. | This routine is provided to clear them (you can also write to _fpCCR,
  205. | since it is globally visible).
  206.  
  207.     .globl  SYM (__clear_sticky_bit)
  208.  
  209.     .text
  210.     .even
  211.  
  212. | void __clear_sticky_bits(void);
  213. SYM (__clear_sticky_bit):        
  214.     lea    SYM (_fpCCR),a0
  215.     movew    IMM (0),a0@(STICK)
  216.     rts
  217.  
  218. |=============================================================================
  219. |                           $_exception_handler
  220. |=============================================================================
  221.  
  222.     .globl  $_exception_handler
  223.  
  224.     .text
  225.     .even
  226.  
  227. | This is the common exit point if an exception occurs.
  228. | NOTE: it is NOT callable from C!
  229. | It expects the exception type in d7, the format (SINGLE_FLOAT,
  230. | DOUBLE_FLOAT or LONG_FLOAT) in d6, and the last operation code in d5.
  231. | It sets the corresponding exception and sticky bits, and the format. 
  232. | Depending on the format if fills the corresponding slots for the 
  233. | operands which produced the exception (all this information is provided
  234. | so if you write your own exception handlers you have enough information
  235. | to deal with the problem).
  236. | Then checks to see if the corresponding exception is trap-enabled, 
  237. | in which case it pushes the address of _fpCCR and traps through 
  238. | trap FPTRAP (15 for the moment).
  239.  
  240. FPTRAP = 15
  241.  
  242. $_exception_handler:
  243.     lea    SYM (_fpCCR),a0
  244.     movew    d7,a0@(EBITS)    | set __exception_bits
  245.     orw    d7,a0@(STICK)    | and __sticky_bits
  246.     movew    d6,a0@(FORMT)    | and __format
  247.     movew    d5,a0@(LASTO)    | and __last_operation
  248.  
  249. | Now put the operands in place:
  250.     cmpw    IMM (SINGLE_FLOAT),d6
  251.     beq    1f
  252.     movel    a6@(8),a0@(OPER1)
  253.     movel    a6@(12),a0@(OPER1+4)
  254.     movel    a6@(16),a0@(OPER2)
  255.     movel    a6@(20),a0@(OPER2+4)
  256.     bra    2f
  257. 1:    movel    a6@(8),a0@(OPER1)
  258.     movel    a6@(12),a0@(OPER2)
  259. 2:
  260. | And check whether the exception is trap-enabled:
  261.     andw    a0@(TRAPE),d7    | is exception trap-enabled?
  262.     beq    1f        | no, exit
  263.     pea    SYM (_fpCCR)    | yes, push address of _fpCCR
  264.     trap    IMM (FPTRAP)    | and trap
  265. 1:    moveml    sp@+,d2-d7    | restore data registers
  266.     unlk    a6        | and return
  267.     rts
  268. #endif /* L_floatex */
  269.  
  270. #ifdef  L_mulsi3
  271.     .text
  272.     .proc
  273.     .globl    SYM (__mulsi3)
  274. SYM (__mulsi3):
  275.     movew    sp@(4), d0    /* x0 -> d0 */
  276.     muluw    sp@(10), d0    /* x0*y1 */
  277.     movew    sp@(6), d1    /* x1 -> d1 */
  278.     muluw    sp@(8), d1    /* x1*y0 */
  279.     addw    d1, d0
  280.     swap    d0
  281.     clrw    d0
  282.     movew    sp@(6), d1    /* x1 -> d1 */
  283.     muluw    sp@(10), d1    /* x1*y1 */
  284.     addl    d1, d0
  285.  
  286.     rts
  287. #endif /* L_mulsi3 */
  288.  
  289. #ifdef  L_udivsi3
  290.     .text
  291.     .proc
  292.     .globl    SYM (__udivsi3)
  293. SYM (__udivsi3):
  294.     movel    d2, sp@-
  295.     movel    sp@(12), d1    /* d1 = divisor */
  296.     movel    sp@(8), d0    /* d0 = dividend */
  297.  
  298.     cmpl    IMM (0x10000), d1 /* divisor >= 2 ^ 16 ?   */
  299.     jcc    L3        /* then try next algorithm */
  300.     movel    d0, d2
  301.     clrw    d2
  302.     swap    d2
  303.     divu    d1, d2          /* high quotient in lower word */
  304.     movew    d2, d0        /* save high quotient */
  305.     swap    d0
  306.     movew    sp@(10), d2    /* get low dividend + high rest */
  307.     divu    d1, d2        /* low quotient */
  308.     movew    d2, d0
  309.     jra    L6
  310.  
  311. L3:    movel    d1, d2        /* use d2 as divisor backup */
  312. L4:    lsrl    IMM (1), d1    /* shift divisor */
  313.     lsrl    IMM (1), d0    /* shift dividend */
  314.     cmpl    IMM (0x10000), d1 /* still divisor >= 2 ^ 16 ?  */
  315.     jcc    L4
  316.     divu    d1, d0        /* now we have 16 bit divisor */
  317.     andl    IMM (0xffff), d0 /* mask out divisor, ignore remainder */
  318.  
  319. /* Multiply the 16 bit tentative quotient with the 32 bit divisor.  Because of
  320.    the operand ranges, this might give a 33 bit product.  If this product is
  321.    greater than the dividend, the tentative quotient was too large. */
  322.     movel    d2, d1
  323.     mulu    d0, d1        /* low part, 32 bits */
  324.     swap    d2
  325.     mulu    d0, d2        /* high part, at most 17 bits */
  326.     swap    d2        /* align high part with low part */
  327.     btst    IMM (0), d2    /* high part 17 bits? */
  328.     jne    L5        /* if 17 bits, quotient was too large */
  329.     addl    d2, d1        /* add parts */
  330.     jcs    L5        /* if sum is 33 bits, quotient was too large */
  331.     cmpl    sp@(8), d1    /* compare the sum with the dividend */
  332.     jls    L6        /* if sum > dividend, quotient was too large */
  333. L5:    subql    IMM (1), d0    /* adjust quotient */
  334.  
  335. L6:    movel    sp@+, d2
  336.     rts
  337. #endif /* L_udivsi3 */
  338.  
  339. #ifdef  L_divsi3
  340.     .text
  341.     .proc
  342.     .globl    SYM (__divsi3)
  343. SYM (__divsi3):
  344.     movel    d2, sp@-
  345.  
  346.     moveb    IMM (1), d2    /* sign of result stored in d2 (=1 or =-1) */
  347.     movel    sp@(12), d1    /* d1 = divisor */
  348.     jpl    L1
  349.     negl    d1
  350.     negb    d2        /* change sign because divisor <0  */
  351. L1:    movel    sp@(8), d0    /* d0 = dividend */
  352.     jpl    L2
  353.     negl    d0
  354.     negb    d2
  355.  
  356. L2:    movel    d1, sp@-
  357.     movel    d0, sp@-
  358.     jbsr    SYM (__udivsi3)    /* divide abs(dividend) by abs(divisor) */
  359.     addql    IMM (8), sp
  360.  
  361.     tstb    d2
  362.     jpl    L3
  363.     negl    d0
  364.  
  365. L3:    movel    sp@+, d2
  366.     rts
  367. #endif /* L_divsi3 */
  368.  
  369. #ifdef  L_umodsi3
  370.     .text
  371.     .proc
  372.     .globl    SYM (__umodsi3)
  373. SYM (__umodsi3):
  374.     movel    sp@(8), d1    /* d1 = divisor */
  375.     movel    sp@(4), d0    /* d0 = dividend */
  376.     movel    d1, sp@-
  377.     movel    d0, sp@-
  378.     jbsr    SYM (__udivsi3)
  379.     addql    IMM (8), sp
  380.     movel    sp@(8), d1    /* d1 = divisor */
  381.     movel    d1, sp@-
  382.     movel    d0, sp@-
  383.     jbsr    SYM (__mulsi3)    /* d0 = (a/b)*b */
  384.     addql    IMM (8), sp
  385.     movel    sp@(4), d1    /* d1 = dividend */
  386.     subl    d0, d1        /* d1 = a - (a/b)*b */
  387.     movel    d1, d0
  388.     rts
  389. #endif /* L_umodsi3 */
  390.  
  391. #ifdef  L_modsi3
  392.     .text
  393.     .proc
  394.     .globl    SYM (__modsi3)
  395. SYM (__modsi3):
  396.     movel    sp@(8), d1    /* d1 = divisor */
  397.     movel    sp@(4), d0    /* d0 = dividend */
  398.     movel    d1, sp@-
  399.     movel    d0, sp@-
  400.     jbsr    SYM (__divsi3)
  401.     addql    IMM (8), sp
  402.     movel    sp@(8), d1    /* d1 = divisor */
  403.     movel    d1, sp@-
  404.     movel    d0, sp@-
  405.     jbsr    SYM (__mulsi3)    /* d0 = (a/b)*b */
  406.     addql    IMM (8), sp
  407.     movel    sp@(4), d1    /* d1 = dividend */
  408.     subl    d0, d1        /* d1 = a - (a/b)*b */
  409.     movel    d1, d0
  410.     rts
  411. #endif /* L_modsi3 */
  412.  
  413.  
  414. #ifdef  L_double
  415.  
  416.     .globl    SYM (_fpCCR)
  417.     .globl  $_exception_handler
  418.  
  419. QUIET_NaN      = 0xffffffff
  420.  
  421. D_MAX_EXP      = 0x07ff
  422. D_BIAS         = 1022
  423. DBL_MAX_EXP    = D_MAX_EXP - D_BIAS
  424. DBL_MIN_EXP    = 1 - D_BIAS
  425. DBL_MANT_DIG   = 53
  426.  
  427. INEXACT_RESULT         = 0x0001
  428. UNDERFLOW         = 0x0002
  429. OVERFLOW         = 0x0004
  430. DIVIDE_BY_ZERO         = 0x0008
  431. INVALID_OPERATION     = 0x0010
  432.  
  433. DOUBLE_FLOAT = 2
  434.  
  435. NOOP         = 0
  436. ADD          = 1
  437. MULTIPLY     = 2
  438. DIVIDE       = 3
  439. NEGATE       = 4
  440. COMPARE      = 5
  441. EXTENDSFDF   = 6
  442. TRUNCDFSF    = 7
  443.  
  444. UNKNOWN           = -1
  445. ROUND_TO_NEAREST  = 0 | round result to nearest representable value
  446. ROUND_TO_ZERO     = 1 | round result towards zero
  447. ROUND_TO_PLUS     = 2 | round result towards plus infinity
  448. ROUND_TO_MINUS    = 3 | round result towards minus infinity
  449.  
  450. | Entry points:
  451.  
  452.     .globl SYM (__adddf3)
  453.     .globl SYM (__subdf3)
  454.     .globl SYM (__muldf3)
  455.     .globl SYM (__divdf3)
  456.     .globl SYM (__negdf2)
  457.     .globl SYM (__cmpdf2)
  458.  
  459.     .text
  460.     .even
  461.  
  462. | These are common routines to return and signal exceptions.    
  463.  
  464. Ld$den:
  465. | Return and signal a denormalized number
  466.     orl    d7,d0
  467.     movew    IMM (UNDERFLOW),d7
  468.     orw    IMM (INEXACT_RESULT),d7
  469.     movew    IMM (DOUBLE_FLOAT),d6
  470.     jmp    $_exception_handler
  471.  
  472. Ld$infty:
  473. Ld$overflow:
  474. | Return a properly signed INFINITY and set the exception flags 
  475.     movel    IMM (0x7ff00000),d0
  476.     movel    IMM (0),d1
  477.     orl    d7,d0
  478.     movew    IMM (OVERFLOW),d7
  479.     orw    IMM (INEXACT_RESULT),d7
  480.     movew    IMM (DOUBLE_FLOAT),d6
  481.     jmp    $_exception_handler
  482.  
  483. Ld$underflow:
  484. | Return 0 and set the exception flags 
  485.     movel    IMM (0),d0
  486.     movel    d0,d1
  487.     movew    IMM (UNDERFLOW),d7
  488.     orw    IMM (INEXACT_RESULT),d7
  489.     movew    IMM (DOUBLE_FLOAT),d6
  490.     jmp    $_exception_handler
  491.  
  492. Ld$inop:
  493. | Return a quiet NaN and set the exception flags
  494.     movel    IMM (QUIET_NaN),d0
  495.     movel    d0,d1
  496.     movew    IMM (INVALID_OPERATION),d7
  497.     orw    IMM (INEXACT_RESULT),d7
  498.     movew    IMM (DOUBLE_FLOAT),d6
  499.     jmp    $_exception_handler
  500.  
  501. Ld$div$0:
  502. | Return a properly signed INFINITY and set the exception flags
  503.     movel    IMM (0x7ff00000),d0
  504.     movel    IMM (0),d1
  505.     orl    d7,d0
  506.     movew    IMM (DIVIDE_BY_ZERO),d7
  507.     orw    IMM (INEXACT_RESULT),d7
  508.     movew    IMM (DOUBLE_FLOAT),d6
  509.     jmp    $_exception_handler
  510.  
  511. |=============================================================================
  512. |=============================================================================
  513. |                         double precision routines
  514. |=============================================================================
  515. |=============================================================================
  516.  
  517. | A double precision floating point number (double) has the format:
  518. |
  519. | struct _double {
  520. |  unsigned int sign      : 1;  /* sign bit */ 
  521. |  unsigned int exponent  : 11; /* exponent, shifted by 126 */
  522. |  unsigned int fraction  : 52; /* fraction */
  523. | } double;
  524. | 
  525. | Thus sizeof(double) = 8 (64 bits). 
  526. |
  527. | All the routines are callable from C programs, and return the result 
  528. | in the register pair d0-d1. They also preserve all registers except 
  529. | d0-d1 and a0-a1.
  530.  
  531. |=============================================================================
  532. |                              __subdf3
  533. |=============================================================================
  534.  
  535. | double __subdf3(double, double);
  536. SYM (__subdf3):
  537.     bchg    IMM (31),sp@(12) | change sign of second operand
  538.                 | and fall through, so we always add
  539. |=============================================================================
  540. |                              __adddf3
  541. |=============================================================================
  542.  
  543. | double __adddf3(double, double);
  544. SYM (__adddf3):
  545.     link    a6,IMM (0)    | everything will be done in registers
  546.     moveml    d2-d7,sp@-    | save all data registers and a2 (but d0-d1)
  547.     movel    a6@(8),d0    | get first operand
  548.     movel    a6@(12),d1    | 
  549.     movel    a6@(16),d2    | get second operand
  550.     movel    a6@(20),d3    | 
  551.  
  552.     movel    d0,d7        | get d0's sign bit in d7 '
  553.     addl    d1,d1        | check and clear sign bit of a, and gain one
  554.     addxl    d0,d0        | bit of extra precision
  555.     beq    Ladddf$b    | if zero return second operand
  556.  
  557.     movel    d2,d6        | save sign in d6 
  558.     addl    d3,d3        | get rid of sign bit and gain one bit of
  559.     addxl    d2,d2        | extra precision
  560.     beq    Ladddf$a    | if zero return first operand
  561.  
  562.     andl    IMM (0x80000000),d7 | isolate a's sign bit '
  563.         swap    d6        | and also b's sign bit '
  564.     andw    IMM (0x8000),d6    |
  565.     orw    d6,d7        | and combine them into d7, so that a's sign '
  566.                 | bit is in the high word and b's is in the '
  567.                 | low word, so d6 is free to be used
  568.     movel    d7,a0        | now save d7 into a0, so d7 is free to
  569.                         | be used also
  570.  
  571. | Get the exponents and check for denormalized and/or infinity.
  572.  
  573.     movel    IMM (0x001fffff),d6 | mask for the fraction
  574.     movel    IMM (0x00200000),d7 | mask to put hidden bit back
  575.  
  576.     movel    d0,d4        | 
  577.     andl    d6,d0        | get fraction in d0
  578.     notl    d6        | make d6 into mask for the exponent
  579.     andl    d6,d4        | get exponent in d4
  580.     beq    Ladddf$a$den    | branch if a is denormalized
  581.     cmpl    d6,d4        | check for INFINITY or NaN
  582.     beq    Ladddf$nf       | 
  583.     orl    d7,d0        | and put hidden bit back
  584. Ladddf$1:
  585.     swap    d4        | shift right exponent so that it starts
  586.     lsrw    IMM (5),d4    | in bit 0 and not bit 20
  587. | Now we have a's exponent in d4 and fraction in d0-d1 '
  588.     movel    d2,d5        | save b to get exponent
  589.     andl    d6,d5        | get exponent in d5
  590.     beq    Ladddf$b$den    | branch if b is denormalized
  591.     cmpl    d6,d5        | check for INFINITY or NaN
  592.     beq    Ladddf$nf
  593.     notl    d6        | make d6 into mask for the fraction again
  594.     andl    d6,d2        | and get fraction in d2
  595.     orl    d7,d2        | and put hidden bit back
  596. Ladddf$2:
  597.     swap    d5        | shift right exponent so that it starts
  598.     lsrw    IMM (5),d5    | in bit 0 and not bit 20
  599.  
  600. | Now we have b's exponent in d5 and fraction in d2-d3. '
  601.  
  602. | The situation now is as follows: the signs are combined in a0, the 
  603. | numbers are in d0-d1 (a) and d2-d3 (b), and the exponents in d4 (a)
  604. | and d5 (b). To do the rounding correctly we need to keep all the
  605. | bits until the end, so we need to use d0-d1-d2-d3 for the first number
  606. | and d4-d5-d6-d7 for the second. To do this we store (temporarily) the
  607. | exponents in a2-a3.
  608.  
  609.     moveml    a2-a3,sp@-    | save the address registers
  610.  
  611.     movel    d4,a2        | save the exponents
  612.     movel    d5,a3        | 
  613.  
  614.     movel    IMM (0),d7    | and move the numbers around
  615.     movel    d7,d6        |
  616.     movel    d3,d5        |
  617.     movel    d2,d4        |
  618.     movel    d7,d3        |
  619.     movel    d7,d2        |
  620.  
  621. | Here we shift the numbers until the exponents are the same, and put 
  622. | the largest exponent in a2.
  623.     exg    d4,a2        | get exponents back
  624.     exg    d5,a3        |
  625.     cmpw    d4,d5        | compare the exponents
  626.     beq    Ladddf$3    | if equal don't shift '
  627.     bhi    9f        | branch if second exponent is higher
  628.  
  629. | Here we have a's exponent larger than b's, so we have to shift b. We do 
  630. | this by using as counter d2:
  631. 1:    movew    d4,d2        | move largest exponent to d2
  632.     subw    d5,d2        | and subtract second exponent
  633.     exg    d4,a2        | get back the longs we saved
  634.     exg    d5,a3        |
  635. | if difference is too large we don't shift (actually, we can just exit) '
  636.     cmpw    IMM (DBL_MANT_DIG+2),d2
  637.     bge    Ladddf$b$small
  638.     cmpw    IMM (32),d2    | if difference >= 32, shift by longs
  639.     bge    5f
  640. 2:    cmpw    IMM (16),d2    | if difference >= 16, shift by words    
  641.     bge    6f
  642.     bra    3f        | enter dbra loop
  643.  
  644. 4:    lsrl    IMM (1),d4
  645.     roxrl    IMM (1),d5
  646.     roxrl    IMM (1),d6
  647.     roxrl    IMM (1),d7
  648. 3:    dbra    d2,4b
  649.     movel    IMM (0),d2
  650.     movel    d2,d3    
  651.     bra    Ladddf$4
  652. 5:
  653.     movel    d6,d7
  654.     movel    d5,d6
  655.     movel    d4,d5
  656.     movel    IMM (0),d4
  657.     subw    IMM (32),d2
  658.     bra    2b
  659. 6:
  660.     movew    d6,d7
  661.     swap    d7
  662.     movew    d5,d6
  663.     swap    d6
  664.     movew    d4,d5
  665.     swap    d5
  666.     movew    IMM (0),d4
  667.     swap    d4
  668.     subw    IMM (16),d2
  669.     bra    3b
  670.     
  671. 9:    exg    d4,d5
  672.     movew    d4,d6
  673.     subw    d5,d6        | keep d5 (largest exponent) in d4
  674.     exg    d4,a2
  675.     exg    d5,a3
  676. | if difference is too large we don't shift (actually, we can just exit) '
  677.     cmpw    IMM (DBL_MANT_DIG+2),d6
  678.     bge    Ladddf$a$small
  679.     cmpw    IMM (32),d6    | if difference >= 32, shift by longs
  680.     bge    5f
  681. 2:    cmpw    IMM (16),d6    | if difference >= 16, shift by words    
  682.     bge    6f
  683.     bra    3f        | enter dbra loop
  684.  
  685. 4:    lsrl    IMM (1),d0
  686.     roxrl    IMM (1),d1
  687.     roxrl    IMM (1),d2
  688.     roxrl    IMM (1),d3
  689. 3:    dbra    d6,4b
  690.     movel    IMM (0),d7
  691.     movel    d7,d6
  692.     bra    Ladddf$4
  693. 5:
  694.     movel    d2,d3
  695.     movel    d1,d2
  696.     movel    d0,d1
  697.     movel    IMM (0),d0
  698.     subw    IMM (32),d6
  699.     bra    2b
  700. 6:
  701.     movew    d2,d3
  702.     swap    d3
  703.     movew    d1,d2
  704.     swap    d2
  705.     movew    d0,d1
  706.     swap    d1
  707.     movew    IMM (0),d0
  708.     swap    d0
  709.     subw    IMM (16),d6
  710.     bra    3b
  711. Ladddf$3:
  712.     exg    d4,a2    
  713.     exg    d5,a3
  714. Ladddf$4:    
  715. | Now we have the numbers in d0--d3 and d4--d7, the exponent in a2, and
  716. | the signs in a4.
  717.  
  718. | Here we have to decide whether to add or subtract the numbers:
  719.     exg    d7,a0        | get the signs 
  720.     exg    d6,a3        | a3 is free to be used
  721.     movel    d7,d6        |
  722.     movew    IMM (0),d7    | get a's sign in d7 '
  723.     swap    d6              |
  724.     movew    IMM (0),d6    | and b's sign in d6 '
  725.     eorl    d7,d6        | compare the signs
  726.     bmi    Lsubdf$0    | if the signs are different we have 
  727.                 | to subtract
  728.     exg    d7,a0        | else we add the numbers
  729.     exg    d6,a3        |
  730.     addl    d7,d3        |
  731.     addxl    d6,d2        |
  732.     addxl    d5,d1        | 
  733.     addxl    d4,d0           |
  734.  
  735.     movel    a2,d4        | return exponent to d4
  736.     movel    a0,d7        | 
  737.     andl    IMM (0x80000000),d7 | d7 now has the sign
  738.  
  739.     moveml    sp@+,a2-a3    
  740.  
  741. | Before rounding normalize so bit #DBL_MANT_DIG is set (we will consider
  742. | the case of denormalized numbers in the rounding routine itself).
  743. | As in the addition (not in the subtraction!) we could have set 
  744. | one more bit we check this:
  745.     btst    IMM (DBL_MANT_DIG+1),d0    
  746.     beq    1f
  747.     lsrl    IMM (1),d0
  748.     roxrl    IMM (1),d1
  749.     roxrl    IMM (1),d2
  750.     roxrl    IMM (1),d3
  751.     addw    IMM (1),d4
  752. 1:
  753.     lea    Ladddf$5,a0    | to return from rounding routine
  754.     lea    SYM (_fpCCR),a1    | check the rounding mode
  755.     movew    a1@(6),d6    | rounding mode in d6
  756.     beq    Lround$to$nearest
  757.     cmpw    IMM (ROUND_TO_PLUS),d6
  758.     bhi    Lround$to$minus
  759.     blt    Lround$to$zero
  760.     bra    Lround$to$plus
  761. Ladddf$5:
  762. | Put back the exponent and check for overflow
  763.     cmpw    IMM (0x7ff),d4    | is the exponent big?
  764.     bge    1f
  765.     bclr    IMM (DBL_MANT_DIG-1),d0
  766.     lslw    IMM (4),d4    | put exponent back into position
  767.     swap    d0        | 
  768.     orw    d4,d0        |
  769.     swap    d0        |
  770.     bra    Ladddf$ret
  771. 1:
  772.     movew    IMM (ADD),d5
  773.     bra    Ld$overflow
  774.  
  775. Lsubdf$0:
  776. | Here we do the subtraction.
  777.     exg    d7,a0        | put sign back in a0
  778.     exg    d6,a3        |
  779.     subl    d7,d3        |
  780.     subxl    d6,d2        |
  781.     subxl    d5,d1        |
  782.     subxl    d4,d0        |
  783.     beq    Ladddf$ret$1    | if zero just exit
  784.     bpl    1f        | if positive skip the following
  785.     exg    d7,a0        |
  786.     bchg    IMM (31),d7    | change sign bit in d7
  787.     exg    d7,a0        |
  788.     negl    d3        |
  789.     negxl    d2        |
  790.     negxl    d1              | and negate result
  791.     negxl    d0              |
  792. 1:    
  793.     movel    a2,d4        | return exponent to d4
  794.     movel    a0,d7
  795.     andl    IMM (0x80000000),d7 | isolate sign bit
  796.     moveml    sp@+,a2-a3    |
  797.  
  798. | Before rounding normalize so bit #DBL_MANT_DIG is set (we will consider
  799. | the case of denormalized numbers in the rounding routine itself).
  800. | As in the addition (not in the subtraction!) we could have set 
  801. | one more bit we check this:
  802.     btst    IMM (DBL_MANT_DIG+1),d0    
  803.     beq    1f
  804.     lsrl    IMM (1),d0
  805.     roxrl    IMM (1),d1
  806.     roxrl    IMM (1),d2
  807.     roxrl    IMM (1),d3
  808.     addw    IMM (1),d4
  809. 1:
  810.     lea    Lsubdf$1,a0    | to return from rounding routine
  811.     lea    SYM (_fpCCR),a1    | check the rounding mode
  812.     movew    a1@(6),d6    | rounding mode in d6
  813.     beq    Lround$to$nearest
  814.     cmpw    IMM (ROUND_TO_PLUS),d6
  815.     bhi    Lround$to$minus
  816.     blt    Lround$to$zero
  817.     bra    Lround$to$plus
  818. Lsubdf$1:
  819. | Put back the exponent and sign (we don't have overflow). '
  820.     bclr    IMM (DBL_MANT_DIG-1),d0    
  821.     lslw    IMM (4),d4    | put exponent back into position
  822.     swap    d0        | 
  823.     orw    d4,d0        |
  824.     swap    d0        |
  825.     bra    Ladddf$ret
  826.  
  827. | If one of the numbers was too small (difference of exponents >= 
  828. | DBL_MANT_DIG+1) we return the other (and now we don't have to '
  829. | check for finiteness or zero).
  830. Ladddf$a$small:
  831.     moveml    sp@+,a2-a3    
  832.     movel    a6@(16),d0
  833.     movel    a6@(20),d1
  834.     lea    SYM (_fpCCR),a0
  835.     movew    IMM (0),a0@
  836.     moveml    sp@+,d2-d7    | restore data registers
  837.     unlk    a6        | and return
  838.     rts
  839.  
  840. Ladddf$b$small:
  841.     moveml    sp@+,a2-a3    
  842.     movel    a6@(8),d0
  843.     movel    a6@(12),d1
  844.     lea    SYM (_fpCCR),a0
  845.     movew    IMM (0),a0@
  846.     moveml    sp@+,d2-d7    | restore data registers
  847.     unlk    a6        | and return
  848.     rts
  849.  
  850. Ladddf$a$den:
  851.     movel    d7,d4        | d7 contains 0x00200000
  852.     bra    Ladddf$1
  853.  
  854. Ladddf$b$den:
  855.     movel    d7,d5           | d7 contains 0x00200000
  856.     notl    d6
  857.     bra    Ladddf$2
  858.  
  859. Ladddf$b:
  860. | Return b (if a is zero)
  861.     movel    d2,d0
  862.     movel    d3,d1
  863.     bra    1f
  864. Ladddf$a:
  865.     movel    a6@(8),d0
  866.     movel    a6@(12),d1
  867. 1:
  868.     movew    IMM (ADD),d5
  869. | Check for NaN and +/-INFINITY.
  870.     movel    d0,d7                 |
  871.     andl    IMM (0x80000000),d7    |
  872.     bclr    IMM (31),d0        |
  873.     cmpl    IMM (0x7ff00000),d0    |
  874.     bge    2f            |
  875.     movel    d0,d0               | check for zero, since we don't  '
  876.     bne    Ladddf$ret        | want to return -0 by mistake
  877.     bclr    IMM (31),d7        |
  878.     bra    Ladddf$ret        |
  879. 2:
  880.     andl    IMM (0x000fffff),d0    | check for NaN (nonzero fraction)
  881.     orl    d1,d0            |
  882.     bne    Ld$inop             |
  883.     bra    Ld$infty        |
  884.     
  885. Ladddf$ret$1:
  886.     moveml    sp@+,a2-a3    | restore regs and exit
  887.  
  888. Ladddf$ret:
  889. | Normal exit.
  890.     lea    SYM (_fpCCR),a0
  891.     movew    IMM (0),a0@
  892.     orl    d7,d0        | put sign bit back
  893.     moveml    sp@+,d2-d7
  894.     unlk    a6
  895.     rts
  896.  
  897. Ladddf$ret$den:
  898. | Return a denormalized number.
  899.     lsrl    IMM (1),d0    | shift right once more
  900.     roxrl    IMM (1),d1    |
  901.     bra    Ladddf$ret
  902.  
  903. Ladddf$nf:
  904.     movew    IMM (ADD),d5
  905. | This could be faster but it is not worth the effort, since it is not
  906. | executed very often. We sacrifice speed for clarity here.
  907.     movel    a6@(8),d0    | get the numbers back (remember that we
  908.     movel    a6@(12),d1    | did some processing already)
  909.     movel    a6@(16),d2    | 
  910.     movel    a6@(20),d3    | 
  911.     movel    IMM (0x7ff00000),d4 | useful constant (INFINITY)
  912.     movel    d0,d7        | save sign bits
  913.     movel    d2,d6        | 
  914.     bclr    IMM (31),d0    | clear sign bits
  915.     bclr    IMM (31),d2    | 
  916. | We know that one of them is either NaN of +/-INFINITY
  917. | Check for NaN (if either one is NaN return NaN)
  918.     cmpl    d4,d0        | check first a (d0)
  919.     bhi    Ld$inop        | if d0 > 0x7ff00000 or equal and
  920.     bne    2f
  921.     tstl    d1        | d1 > 0, a is NaN
  922.     bne    Ld$inop        | 
  923. 2:    cmpl    d4,d2        | check now b (d1)
  924.     bhi    Ld$inop        | 
  925.     bne    3f
  926.     tstl    d3        | 
  927.     bne    Ld$inop        | 
  928. 3:
  929. | Now comes the check for +/-INFINITY. We know that both are (maybe not
  930. | finite) numbers, but we have to check if both are infinite whether we
  931. | are adding or subtracting them.
  932.     eorl    d7,d6        | to check sign bits
  933.     bmi    1f
  934.     andl    IMM (0x80000000),d7 | get (common) sign bit
  935.     bra    Ld$infty
  936. 1:
  937. | We know one (or both) are infinite, so we test for equality between the
  938. | two numbers (if they are equal they have to be infinite both, so we
  939. | return NaN).
  940.     cmpl    d2,d0        | are both infinite?
  941.     bne    1f        | if d0 <> d2 they are not equal
  942.     cmpl    d3,d1        | if d0 == d2 test d3 and d1
  943.     beq    Ld$inop        | if equal return NaN
  944. 1:    
  945.     andl    IMM (0x80000000),d7 | get a's sign bit '
  946.     cmpl    d4,d0        | test now for infinity
  947.     beq    Ld$infty    | if a is INFINITY return with this sign
  948.     bchg    IMM (31),d7    | else we know b is INFINITY and has
  949.     bra    Ld$infty    | the opposite sign
  950.  
  951. |=============================================================================
  952. |                              __muldf3
  953. |=============================================================================
  954.  
  955. | double __muldf3(double, double);
  956. SYM (__muldf3):
  957.     link    a6,IMM (0)
  958.     moveml    d2-d7,sp@-
  959.     movel    a6@(8),d0        | get a into d0-d1
  960.     movel    a6@(12),d1        | 
  961.     movel    a6@(16),d2        | and b into d2-d3
  962.     movel    a6@(20),d3        |
  963.     movel    d0,d7            | d7 will hold the sign of the product
  964.     eorl    d2,d7            |
  965.     andl    IMM (0x80000000),d7    |
  966.     movel    d7,a0            | save sign bit into a0 
  967.     movel    IMM (0x7ff00000),d7    | useful constant (+INFINITY)
  968.     movel    d7,d6            | another (mask for fraction)
  969.     notl    d6            |
  970.     bclr    IMM (31),d0        | get rid of a's sign bit '
  971.     movel    d0,d4            | 
  972.     orl    d1,d4            | 
  973.     beq    Lmuldf$a$0        | branch if a is zero
  974.     movel    d0,d4            |
  975.     bclr    IMM (31),d2        | get rid of b's sign bit '
  976.     movel    d2,d5            |
  977.     orl    d3,d5            | 
  978.     beq    Lmuldf$b$0        | branch if b is zero
  979.     movel    d2,d5            | 
  980.     cmpl    d7,d0            | is a big?
  981.     bhi    Lmuldf$inop        | if a is NaN return NaN
  982.     beq    Lmuldf$a$nf        | we still have to check d1 and b ...
  983.     cmpl    d7,d2            | now compare b with INFINITY
  984.     bhi    Lmuldf$inop        | is b NaN?
  985.     beq    Lmuldf$b$nf         | we still have to check d3 ...
  986. | Here we have both numbers finite and nonzero (and with no sign bit).
  987. | Now we get the exponents into d4 and d5.
  988.     andl    d7,d4            | isolate exponent in d4
  989.     beq    Lmuldf$a$den        | if exponent zero, have denormalized
  990.     andl    d6,d0            | isolate fraction
  991.     orl    IMM (0x00100000),d0    | and put hidden bit back
  992.     swap    d4            | I like exponents in the first byte
  993.     lsrw    IMM (4),d4        | 
  994. Lmuldf$1:            
  995.     andl    d7,d5            |
  996.     beq    Lmuldf$b$den        |
  997.     andl    d6,d2            |
  998.     orl    IMM (0x00100000),d2    | and put hidden bit back
  999.     swap    d5            |
  1000.     lsrw    IMM (4),d5        |
  1001. Lmuldf$2:                |
  1002.     addw    d5,d4            | add exponents
  1003.     subw    IMM (D_BIAS+1),d4    | and subtract bias (plus one)
  1004.  
  1005. | We are now ready to do the multiplication. The situation is as follows:
  1006. | both a and b have bit 52 ( bit 20 of d0 and d2) set (even if they were 
  1007. | denormalized to start with!), which means that in the product bit 104 
  1008. | (which will correspond to bit 8 of the fourth long) is set.
  1009.  
  1010. | Here we have to do the product.
  1011. | To do it we have to juggle the registers back and forth, as there are not
  1012. | enough to keep everything in them. So we use the address registers to keep
  1013. | some intermediate data.
  1014.  
  1015.     moveml    a2-a3,sp@-    | save a2 and a3 for temporary use
  1016.     movel    IMM (0),a2    | a2 is a null register
  1017.     movel    d4,a3        | and a3 will preserve the exponent
  1018.  
  1019. | First, shift d2-d3 so bit 20 becomes bit 31:
  1020.     rorl    IMM (5),d2    | rotate d2 5 places right
  1021.     swap    d2        | and swap it
  1022.     rorl    IMM (5),d3    | do the same thing with d3
  1023.     swap    d3        |
  1024.     movew    d3,d6        | get the rightmost 11 bits of d3
  1025.     andw    IMM (0x07ff),d6    |
  1026.     orw    d6,d2        | and put them into d2
  1027.     andw    IMM (0xf800),d3    | clear those bits in d3
  1028.  
  1029.     movel    d2,d6        | move b into d6-d7
  1030.     movel    d3,d7           | move a into d4-d5
  1031.     movel    d0,d4           | and clear d0-d1-d2-d3 (to put result)
  1032.     movel    d1,d5           |
  1033.     movel    IMM (0),d3    |
  1034.     movel    d3,d2           |
  1035.     movel    d3,d1           |
  1036.     movel    d3,d0            |
  1037.  
  1038. | We use a1 as counter:    
  1039.     movel    IMM (DBL_MANT_DIG-1),a1        
  1040.     exg    d7,a1
  1041.  
  1042. 1:    exg    d7,a1        | put counter back in a1
  1043.     addl    d3,d3        | shift sum once left
  1044.     addxl    d2,d2           |
  1045.     addxl    d1,d1           |
  1046.     addxl    d0,d0           |
  1047.     addl    d7,d7        |
  1048.     addxl    d6,d6        |
  1049.     bcc    2f        | if bit clear skip the following
  1050.     exg    d7,a2        |
  1051.     addl    d5,d3        | else add a to the sum
  1052.     addxl    d4,d2        |
  1053.     addxl    d7,d1        |
  1054.     addxl    d7,d0        |
  1055.     exg    d7,a2        | 
  1056. 2:    exg    d7,a1        | put counter in d7
  1057.     dbf    d7,1b        | decrement and branch
  1058.  
  1059.     movel    a3,d4        | restore exponent
  1060.     moveml    sp@+,a2-a3
  1061.  
  1062. | Now we have the product in d0-d1-d2-d3, with bit 8 of d0 set. The 
  1063. | first thing to do now is to normalize it so bit 8 becomes bit 
  1064. | DBL_MANT_DIG-32 (to do the rounding); later we will shift right.
  1065.     swap    d0
  1066.     swap    d1
  1067.     movew    d1,d0
  1068.     swap    d2
  1069.     movew    d2,d1
  1070.     swap    d3
  1071.     movew    d3,d2
  1072.     movew    IMM (0),d3
  1073.     lsrl    IMM (1),d0
  1074.     roxrl    IMM (1),d1
  1075.     roxrl    IMM (1),d2
  1076.     roxrl    IMM (1),d3
  1077.     lsrl    IMM (1),d0
  1078.     roxrl    IMM (1),d1
  1079.     roxrl    IMM (1),d2
  1080.     roxrl    IMM (1),d3
  1081.     lsrl    IMM (1),d0
  1082.     roxrl    IMM (1),d1
  1083.     roxrl    IMM (1),d2
  1084.     roxrl    IMM (1),d3
  1085.     
  1086. | Now round, check for over- and underflow, and exit.
  1087.     movel    a0,d7        | get sign bit back into d7
  1088.     movew    IMM (MULTIPLY),d5
  1089.  
  1090.     btst    IMM (DBL_MANT_DIG+1-32),d0
  1091.     beq    Lround$exit
  1092.     lsrl    IMM (1),d0
  1093.     roxrl    IMM (1),d1
  1094.     addw    IMM (1),d4
  1095.     bra    Lround$exit
  1096.  
  1097. Lmuldf$inop:
  1098.     movew    IMM (MULTIPLY),d5
  1099.     bra    Ld$inop
  1100.  
  1101. Lmuldf$b$nf:
  1102.     movew    IMM (MULTIPLY),d5
  1103.     movel    a0,d7        | get sign bit back into d7
  1104.     tstl    d3        | we know d2 == 0x7ff00000, so check d3
  1105.     bne    Ld$inop        | if d3 <> 0 b is NaN
  1106.     bra    Ld$overflow    | else we have overflow (since a is finite)
  1107.  
  1108. Lmuldf$a$nf:
  1109.     movew    IMM (MULTIPLY),d5
  1110.     movel    a0,d7        | get sign bit back into d7
  1111.     tstl    d1        | we know d0 == 0x7ff00000, so check d1
  1112.     bne    Ld$inop        | if d1 <> 0 a is NaN
  1113.     bra    Ld$overflow    | else signal overflow
  1114.  
  1115. | If either number is zero return zero, unless the other is +/-INFINITY or
  1116. | NaN, in which case we return NaN.
  1117. Lmuldf$b$0:
  1118.     movew    IMM (MULTIPLY),d5
  1119.     exg    d2,d0        | put b (==0) into d0-d1
  1120.     exg    d3,d1        | and a (with sign bit cleared) into d2-d3
  1121.     bra    1f
  1122. Lmuldf$a$0:
  1123.     movel    a6@(16),d2    | put b into d2-d3 again
  1124.     movel    a6@(20),d3    |
  1125.     bclr    IMM (31),d2    | clear sign bit
  1126. 1:    cmpl    IMM (0x7ff00000),d2 | check for non-finiteness
  1127.     bge    Ld$inop        | in case NaN or +/-INFINITY return NaN
  1128.     lea    SYM (_fpCCR),a0
  1129.     movew    IMM (0),a0@
  1130.     moveml    sp@+,d2-d7
  1131.     unlk    a6
  1132.     rts
  1133.  
  1134. | If a number is denormalized we put an exponent of 1 but do not put the 
  1135. | hidden bit back into the fraction; instead we shift left until bit 21
  1136. | (the hidden bit) is set, adjusting the exponent accordingly. We do this
  1137. | to ensure that the product of the fractions is close to 1.
  1138. Lmuldf$a$den:
  1139.     movel    IMM (1),d4
  1140.     andl    d6,d0
  1141. 1:    addl    d1,d1           | shift a left until bit 20 is set
  1142.     addxl    d0,d0        |
  1143.     subw    IMM (1),d4    | and adjust exponent
  1144.     btst    IMM (20),d0    |
  1145.     bne    Lmuldf$1        |
  1146.     bra    1b
  1147.  
  1148. Lmuldf$b$den:
  1149.     movel    IMM (1),d5
  1150.     andl    d6,d2
  1151. 1:    addl    d3,d3        | shift b left until bit 20 is set
  1152.     addxl    d2,d2        |
  1153.     subw    IMM (1),d5    | and adjust exponent
  1154.     btst    IMM (20),d2    |
  1155.     bne    Lmuldf$2    |
  1156.     bra    1b
  1157.  
  1158.  
  1159. |=============================================================================
  1160. |                              __divdf3
  1161. |=============================================================================
  1162.  
  1163. | double __divdf3(double, double);
  1164. SYM (__divdf3):
  1165.     link    a6,IMM (0)
  1166.     moveml    d2-d7,sp@-
  1167.     movel    a6@(8),d0    | get a into d0-d1
  1168.     movel    a6@(12),d1    | 
  1169.     movel    a6@(16),d2    | and b into d2-d3
  1170.     movel    a6@(20),d3    |
  1171.     movel    d0,d7        | d7 will hold the sign of the result
  1172.     eorl    d2,d7        |
  1173.     andl    IMM (0x80000000),d7
  1174.     movel    d7,a0        | save sign into a0
  1175.     movel    IMM (0x7ff00000),d7 | useful constant (+INFINITY)
  1176.     movel    d7,d6        | another (mask for fraction)
  1177.     notl    d6        |
  1178.     bclr    IMM (31),d0    | get rid of a's sign bit '
  1179.     movel    d0,d4        |
  1180.     orl    d1,d4        |
  1181.     beq    Ldivdf$a$0    | branch if a is zero
  1182.     movel    d0,d4        |
  1183.     bclr    IMM (31),d2    | get rid of b's sign bit '
  1184.     movel    d2,d5        |
  1185.     orl    d3,d5        |
  1186.     beq    Ldivdf$b$0    | branch if b is zero
  1187.     movel    d2,d5
  1188.     cmpl    d7,d0        | is a big?
  1189.     bhi    Ldivdf$inop    | if a is NaN return NaN
  1190.     beq    Ldivdf$a$nf    | if d0 == 0x7ff00000 we check d1
  1191.     cmpl    d7,d2        | now compare b with INFINITY 
  1192.     bhi    Ldivdf$inop    | if b is NaN return NaN
  1193.     beq    Ldivdf$b$nf    | if d2 == 0x7ff00000 we check d3
  1194. | Here we have both numbers finite and nonzero (and with no sign bit).
  1195. | Now we get the exponents into d4 and d5 and normalize the numbers to
  1196. | ensure that the ratio of the fractions is around 1. We do this by
  1197. | making sure that both numbers have bit #DBL_MANT_DIG-32-1 (hidden bit)
  1198. | set, even if they were denormalized to start with.
  1199. | Thus, the result will satisfy: 2 > result > 1/2.
  1200.     andl    d7,d4        | and isolate exponent in d4
  1201.     beq    Ldivdf$a$den    | if exponent is zero we have a denormalized
  1202.     andl    d6,d0        | and isolate fraction
  1203.     orl    IMM (0x00100000),d0 | and put hidden bit back
  1204.     swap    d4        | I like exponents in the first byte
  1205.     lsrw    IMM (4),d4    | 
  1206. Ldivdf$1:            | 
  1207.     andl    d7,d5        |
  1208.     beq    Ldivdf$b$den    |
  1209.     andl    d6,d2        |
  1210.     orl    IMM (0x00100000),d2
  1211.     swap    d5        |
  1212.     lsrw    IMM (4),d5    |
  1213. Ldivdf$2:            |
  1214.     subw    d5,d4        | subtract exponents
  1215.     addw    IMM (D_BIAS),d4    | and add bias
  1216.  
  1217. | We are now ready to do the division. We have prepared things in such a way
  1218. | that the ratio of the fractions will be less than 2 but greater than 1/2.
  1219. | At this point the registers in use are:
  1220. | d0-d1    hold a (first operand, bit DBL_MANT_DIG-32=0, bit 
  1221. | DBL_MANT_DIG-1-32=1)
  1222. | d2-d3    hold b (second operand, bit DBL_MANT_DIG-32=1)
  1223. | d4    holds the difference of the exponents, corrected by the bias
  1224. | a0    holds the sign of the ratio
  1225.  
  1226. | To do the rounding correctly we need to keep information about the
  1227. | nonsignificant bits. One way to do this would be to do the division
  1228. | using four registers; another is to use two registers (as originally
  1229. | I did), but use a sticky bit to preserve information about the 
  1230. | fractional part. Note that we can keep that info in a1, which is not
  1231. | used.
  1232.     movel    IMM (0),d6    | d6-d7 will hold the result
  1233.     movel    d6,d7        | 
  1234.     movel    IMM (0),a1    | and a1 will hold the sticky bit
  1235.  
  1236.     movel    IMM (DBL_MANT_DIG-32+1),d5    
  1237.     
  1238. 1:    cmpl    d0,d2        | is a < b?
  1239.     bhi    3f        | if b > a skip the following
  1240.     beq    4f        | if d0==d2 check d1 and d3
  1241. 2:    subl    d3,d1        | 
  1242.     subxl    d2,d0        | a <-- a - b
  1243.     bset    d5,d6        | set the corresponding bit in d6
  1244. 3:    addl    d1,d1        | shift a by 1
  1245.     addxl    d0,d0        |
  1246.     dbra    d5,1b        | and branch back
  1247.     bra    5f            
  1248. 4:    cmpl    d1,d3        | here d0==d2, so check d1 and d3
  1249.     bhi    3b        | if d1 > d2 skip the subtraction
  1250.     bra    2b        | else go do it
  1251. 5:
  1252. | Here we have to start setting the bits in the second long.
  1253.     movel    IMM (31),d5    | again d5 is counter
  1254.  
  1255. 1:    cmpl    d0,d2        | is a < b?
  1256.     bhi    3f        | if b > a skip the following
  1257.     beq    4f        | if d0==d2 check d1 and d3
  1258. 2:    subl    d3,d1        | 
  1259.     subxl    d2,d0        | a <-- a - b
  1260.     bset    d5,d7        | set the corresponding bit in d7
  1261. 3:    addl    d1,d1        | shift a by 1
  1262.     addxl    d0,d0        |
  1263.     dbra    d5,1b        | and branch back
  1264.     bra    5f            
  1265. 4:    cmpl    d1,d3        | here d0==d2, so check d1 and d3
  1266.     bhi    3b        | if d1 > d2 skip the subtraction
  1267.     bra    2b        | else go do it
  1268. 5:
  1269. | Now go ahead checking until we hit a one, which we store in d2.
  1270.     movel    IMM (DBL_MANT_DIG),d5
  1271. 1:    cmpl    d2,d0        | is a < b?
  1272.     bhi    4f        | if b < a, exit
  1273.     beq    3f        | if d0==d2 check d1 and d3
  1274. 2:    addl    d1,d1        | shift a by 1
  1275.     addxl    d0,d0        |
  1276.     dbra    d5,1b        | and branch back
  1277.     movel    IMM (0),d2    | here no sticky bit was found
  1278.     movel    d2,d3
  1279.     bra    5f            
  1280. 3:    cmpl    d1,d3        | here d0==d2, so check d1 and d3
  1281.     bhi    2b        | if d1 > d2 go back
  1282. 4:
  1283. | Here put the sticky bit in d2-d3 (in the position which actually corresponds
  1284. | to it; if you don't do this the algorithm loses in some cases). '
  1285.     movel    IMM (0),d2
  1286.     movel    d2,d3
  1287.     subw    IMM (DBL_MANT_DIG),d5
  1288.     addw    IMM (63),d5
  1289.     cmpw    IMM (31),d5
  1290.     bhi    2f
  1291. 1:    bset    d5,d3
  1292.     bra    5f
  1293.     subw    IMM (32),d5
  1294. 2:    bset    d5,d2
  1295. 5:
  1296. | Finally we are finished! Move the longs in the address registers to
  1297. | their final destination:
  1298.     movel    d6,d0
  1299.     movel    d7,d1
  1300.     movel    IMM (0),d3
  1301.  
  1302. | Here we have finished the division, with the result in d0-d1-d2-d3, with
  1303. | 2^21 <= d6 < 2^23. Thus bit 23 is not set, but bit 22 could be set.
  1304. | If it is not, then definitely bit 21 is set. Normalize so bit 22 is
  1305. | not set:
  1306.     btst    IMM (DBL_MANT_DIG-32+1),d0
  1307.     beq    1f
  1308.     lsrl    IMM (1),d0
  1309.     roxrl    IMM (1),d1
  1310.     roxrl    IMM (1),d2
  1311.     roxrl    IMM (1),d3
  1312.     addw    IMM (1),d4
  1313. 1:
  1314. | Now round, check for over- and underflow, and exit.
  1315.     movel    a0,d7        | restore sign bit to d7
  1316.     movew    IMM (DIVIDE),d5
  1317.     bra    Lround$exit
  1318.  
  1319. Ldivdf$inop:
  1320.     movew    IMM (DIVIDE),d5
  1321.     bra    Ld$inop
  1322.  
  1323. Ldivdf$a$0:
  1324. | If a is zero check to see whether b is zero also. In that case return
  1325. | NaN; then check if b is NaN, and return NaN also in that case. Else
  1326. | return zero.
  1327.     movew    IMM (DIVIDE),d5
  1328.     bclr    IMM (31),d2    |
  1329.     movel    d2,d4        | 
  1330.     orl    d3,d4        | 
  1331.     beq    Ld$inop        | if b is also zero return NaN
  1332.     cmpl    IMM (0x7ff00000),d2 | check for NaN
  1333.     bhi    Ld$inop        | 
  1334.     blt    1f        |
  1335.     tstl    d3        |
  1336.     bne    Ld$inop        |
  1337. 1:    movel    IMM (0),d0    | else return zero
  1338.     movel    d0,d1        | 
  1339.     lea    SYM (_fpCCR),a0    | clear exception flags
  1340.     movew    IMM (0),a0@    |
  1341.     moveml    sp@+,d2-d7    | 
  1342.     unlk    a6        | 
  1343.     rts            |     
  1344.  
  1345. Ldivdf$b$0:
  1346.     movew    IMM (DIVIDE),d5
  1347. | If we got here a is not zero. Check if a is NaN; in that case return NaN,
  1348. | else return +/-INFINITY. Remember that a is in d0 with the sign bit 
  1349. | cleared already.
  1350.     movel    a0,d7        | put a's sign bit back in d7 '
  1351.     cmpl    IMM (0x7ff00000),d0 | compare d0 with INFINITY
  1352.     bhi    Ld$inop        | if larger it is NaN
  1353.     tstl    d1        | 
  1354.     bne    Ld$inop        | 
  1355.     bra    Ld$div$0    | else signal DIVIDE_BY_ZERO
  1356.  
  1357. Ldivdf$b$nf:
  1358.     movew    IMM (DIVIDE),d5
  1359. | If d2 == 0x7ff00000 we have to check d3.
  1360.     tstl    d3        |
  1361.     bne    Ld$inop        | if d3 <> 0, b is NaN
  1362.     bra    Ld$underflow    | else b is +/-INFINITY, so signal underflow
  1363.  
  1364. Ldivdf$a$nf:
  1365.     movew    IMM (DIVIDE),d5
  1366. | If d0 == 0x7ff00000 we have to check d1.
  1367.     tstl    d1        |
  1368.     bne    Ld$inop        | if d1 <> 0, a is NaN
  1369. | If a is INFINITY we have to check b
  1370.     cmpl    d7,d2        | compare b with INFINITY 
  1371.     bge    Ld$inop        | if b is NaN or INFINITY return NaN
  1372.     tstl    d3        |
  1373.     bne    Ld$inop        | 
  1374.     bra    Ld$overflow    | else return overflow
  1375.  
  1376. | If a number is denormalized we put an exponent of 1 but do not put the 
  1377. | bit back into the fraction.
  1378. Ldivdf$a$den:
  1379.     movel    IMM (1),d4
  1380.     andl    d6,d0
  1381. 1:    addl    d1,d1        | shift a left until bit 20 is set
  1382.     addxl    d0,d0
  1383.     subw    IMM (1),d4    | and adjust exponent
  1384.     btst    IMM (DBL_MANT_DIG-32-1),d0
  1385.     bne    Ldivdf$1
  1386.     bra    1b
  1387.  
  1388. Ldivdf$b$den:
  1389.     movel    IMM (1),d5
  1390.     andl    d6,d2
  1391. 1:    addl    d3,d3        | shift b left until bit 20 is set
  1392.     addxl    d2,d2
  1393.     subw    IMM (1),d5    | and adjust exponent
  1394.     btst    IMM (DBL_MANT_DIG-32-1),d2
  1395.     bne    Ldivdf$2
  1396.     bra    1b
  1397.  
  1398. Lround$exit:
  1399. | This is a common exit point for __muldf3 and __divdf3. When they enter
  1400. | this point the sign of the result is in d7, the result in d0-d1, normalized
  1401. | so that 2^21 <= d0 < 2^22, and the exponent is in the lower byte of d4.
  1402.  
  1403. | First check for underlow in the exponent:
  1404.     cmpw    IMM (-DBL_MANT_DIG-1),d4        
  1405.     blt    Ld$underflow    
  1406. | It could happen that the exponent is less than 1, in which case the 
  1407. | number is denormalized. In this case we shift right and adjust the 
  1408. | exponent until it becomes 1 or the fraction is zero (in the latter case 
  1409. | we signal underflow and return zero).
  1410.     movel    d7,a0        |
  1411.     movel    IMM (0),d6    | use d6-d7 to collect bits flushed right
  1412.     movel    d6,d7        | use d6-d7 to collect bits flushed right
  1413.     cmpw    IMM (1),d4    | if the exponent is less than 1 we 
  1414.     bge    2f        | have to shift right (denormalize)
  1415. 1:    addw    IMM (1),d4    | adjust the exponent
  1416.     lsrl    IMM (1),d0    | shift right once 
  1417.     roxrl    IMM (1),d1    |
  1418.     roxrl    IMM (1),d2    |
  1419.     roxrl    IMM (1),d3    |
  1420.     roxrl    IMM (1),d6    | 
  1421.     roxrl    IMM (1),d7    |
  1422.     cmpw    IMM (1),d4    | is the exponent 1 already?
  1423.     beq    2f        | if not loop back
  1424.     bra    1b              |
  1425.     bra    Ld$underflow    | safety check, shouldn't execute '
  1426. 2:    orl    d6,d2        | this is a trick so we don't lose  '
  1427.     orl    d7,d3        | the bits which were flushed right
  1428.     movel    a0,d7        | get back sign bit into d7
  1429. | Now call the rounding routine (which takes care of denormalized numbers):
  1430.     lea    Lround$0,a0    | to return from rounding routine
  1431.     lea    SYM (_fpCCR),a1    | check the rounding mode
  1432.     movew    a1@(6),d6    | rounding mode in d6
  1433.     beq    Lround$to$nearest
  1434.     cmpw    IMM (ROUND_TO_PLUS),d6
  1435.     bhi    Lround$to$minus
  1436.     blt    Lround$to$zero
  1437.     bra    Lround$to$plus
  1438. Lround$0:
  1439. | Here we have a correctly rounded result (either normalized or denormalized).
  1440.  
  1441. | Here we should have either a normalized number or a denormalized one, and
  1442. | the exponent is necessarily larger or equal to 1 (so we don't have to  '
  1443. | check again for underflow!). We have to check for overflow or for a 
  1444. | denormalized number (which also signals underflow).
  1445. | Check for overflow (i.e., exponent >= 0x7ff).
  1446.     cmpw    IMM (0x07ff),d4
  1447.     bge    Ld$overflow
  1448. | Now check for a denormalized number (exponent==0):
  1449.     movew    d4,d4
  1450.     beq    Ld$den
  1451. 1:
  1452. | Put back the exponents and sign and return.
  1453.     lslw    IMM (4),d4    | exponent back to fourth byte
  1454.     bclr    IMM (DBL_MANT_DIG-32-1),d0
  1455.     swap    d0        | and put back exponent
  1456.     orw    d4,d0        | 
  1457.     swap    d0        |
  1458.     orl    d7,d0        | and sign also
  1459.  
  1460.     lea    SYM (_fpCCR),a0
  1461.     movew    IMM (0),a0@
  1462.     moveml    sp@+,d2-d7
  1463.     unlk    a6
  1464.     rts
  1465.  
  1466. |=============================================================================
  1467. |                              __negdf2
  1468. |=============================================================================
  1469.  
  1470. | double __negdf2(double, double);
  1471. SYM (__negdf2):
  1472.     link    a6,IMM (0)
  1473.     moveml    d2-d7,sp@-
  1474.     movew    IMM (NEGATE),d5
  1475.     movel    a6@(8),d0    | get number to negate in d0-d1
  1476.     movel    a6@(12),d1    |
  1477.     bchg    IMM (31),d0    | negate
  1478.     movel    d0,d2        | make a positive copy (for the tests)
  1479.     bclr    IMM (31),d2    |
  1480.     movel    d2,d4        | check for zero
  1481.     orl    d1,d4        |
  1482.     beq    2f        | if zero (either sign) return +zero
  1483.     cmpl    IMM (0x7ff00000),d2 | compare to +INFINITY
  1484.     blt    1f        | if finite, return
  1485.     bhi    Ld$inop        | if larger (fraction not zero) is NaN
  1486.     tstl    d1        | if d2 == 0x7ff00000 check d1
  1487.     bne    Ld$inop        |
  1488.     movel    d0,d7        | else get sign and return INFINITY
  1489.     andl    IMM (0x80000000),d7
  1490.     bra    Ld$infty        
  1491. 1:    lea    SYM (_fpCCR),a0
  1492.     movew    IMM (0),a0@
  1493.     moveml    sp@+,d2-d7
  1494.     unlk    a6
  1495.     rts
  1496. 2:    bclr    IMM (31),d0
  1497.     bra    1b
  1498.  
  1499. |=============================================================================
  1500. |                              __cmpdf2
  1501. |=============================================================================
  1502.  
  1503. GREATER =  1
  1504. LESS    = -1
  1505. EQUAL   =  0
  1506.  
  1507. | int __cmpdf2(double, double);
  1508. SYM (__cmpdf2):
  1509.     link    a6,IMM (0)
  1510.     moveml    d2-d7,sp@-     | save registers
  1511.     movew    IMM (COMPARE),d5
  1512.     movel    a6@(8),d0    | get first operand
  1513.     movel    a6@(12),d1    |
  1514.     movel    a6@(16),d2    | get second operand
  1515.     movel    a6@(20),d3    |
  1516. | First check if a and/or b are (+/-) zero and in that case clear
  1517. | the sign bit.
  1518.     movel    d0,d6        | copy signs into d6 (a) and d7(b)
  1519.     bclr    IMM (31),d0    | and clear signs in d0 and d2
  1520.     movel    d2,d7        |
  1521.     bclr    IMM (31),d2    |
  1522.     cmpl    IMM (0x7fff0000),d0 | check for a == NaN
  1523.     bhi    Ld$inop        | if d0 > 0x7ff00000, a is NaN
  1524.     beq    Lcmpdf$a$nf    | if equal can be INFINITY, so check d1
  1525.     movel    d0,d4        | copy into d4 to test for zero
  1526.     orl    d1,d4        |
  1527.     beq    Lcmpdf$a$0    |
  1528. Lcmpdf$0:
  1529.     cmpl    IMM (0x7fff0000),d2 | check for b == NaN
  1530.     bhi    Ld$inop        | if d2 > 0x7ff00000, b is NaN
  1531.     beq    Lcmpdf$b$nf    | if equal can be INFINITY, so check d3
  1532.     movel    d2,d4        |
  1533.     orl    d3,d4        |
  1534.     beq    Lcmpdf$b$0    |
  1535. Lcmpdf$1:
  1536. | Check the signs
  1537.     eorl    d6,d7
  1538.     bpl    1f
  1539. | If the signs are not equal check if a >= 0
  1540.     tstl    d6
  1541.     bpl    Lcmpdf$a$gt$b    | if (a >= 0 && b < 0) => a > b
  1542.     bmi    Lcmpdf$b$gt$a    | if (a < 0 && b >= 0) => a < b
  1543. 1:
  1544. | If the signs are equal check for < 0
  1545.     tstl    d6
  1546.     bpl    1f
  1547. | If both are negative exchange them
  1548.     exg    d0,d2
  1549.     exg    d1,d3
  1550. 1:
  1551. | Now that they are positive we just compare them as longs (does this also
  1552. | work for denormalized numbers?).
  1553.     cmpl    d0,d2
  1554.     bhi    Lcmpdf$b$gt$a    | |b| > |a|
  1555.     bne    Lcmpdf$a$gt$b    | |b| < |a|
  1556. | If we got here d0 == d2, so we compare d1 and d3.
  1557.     cmpl    d1,d3
  1558.     bhi    Lcmpdf$b$gt$a    | |b| > |a|
  1559.     bne    Lcmpdf$a$gt$b    | |b| < |a|
  1560. | If we got here a == b.
  1561.     movel    IMM (EQUAL),d0
  1562.     moveml    sp@+,d2-d7     | put back the registers
  1563.     unlk    a6
  1564.     rts
  1565. Lcmpdf$a$gt$b:
  1566.     movel    IMM (GREATER),d0
  1567.     moveml    sp@+,d2-d7     | put back the registers
  1568.     unlk    a6
  1569.     rts
  1570. Lcmpdf$b$gt$a:
  1571.     movel    IMM (LESS),d0
  1572.     moveml    sp@+,d2-d7     | put back the registers
  1573.     unlk    a6
  1574.     rts
  1575.  
  1576. Lcmpdf$a$0:    
  1577.     bclr    IMM (31),d6
  1578.     bra    Lcmpdf$0
  1579. Lcmpdf$b$0:
  1580.     bclr    IMM (31),d7
  1581.     bra    Lcmpdf$1
  1582.  
  1583. Lcmpdf$a$nf:
  1584.     tstl    d1
  1585.     bne    Ld$inop
  1586.     bra    Lcmpdf$0
  1587.  
  1588. Lcmpdf$b$nf:
  1589.     tstl    d3
  1590.     bne    Ld$inop
  1591.     bra    Lcmpdf$1
  1592.  
  1593. |=============================================================================
  1594. |                           rounding routines
  1595. |=============================================================================
  1596.  
  1597. | The rounding routines expect the number to be normalized in registers
  1598. | d0-d1-d2-d3, with the exponent in register d4. They assume that the 
  1599. | exponent is larger or equal to 1. They return a properly normalized number
  1600. | if possible, and a denormalized number otherwise. The exponent is returned
  1601. | in d4.
  1602.  
  1603. Lround$to$nearest:
  1604. | We now normalize as suggested by D. Knuth ("Seminumerical Algorithms"):
  1605. | Here we assume that the exponent is not too small (this should be checked
  1606. | before entering the rounding routine), but the number could be denormalized.
  1607.  
  1608. | Check for denormalized numbers:
  1609. 1:    btst    IMM (DBL_MANT_DIG-32),d0
  1610.     bne    2f        | if set the number is normalized
  1611. | Normalize shifting left until bit #DBL_MANT_DIG-32 is set or the exponent 
  1612. | is one (remember that a denormalized number corresponds to an 
  1613. | exponent of -D_BIAS+1).
  1614.     cmpw    IMM (1),d4    | remember that the exponent is at least one
  1615.      beq    2f        | an exponent of one means denormalized
  1616.     addl    d3,d3        | else shift and adjust the exponent
  1617.     addxl    d2,d2        |
  1618.     addxl    d1,d1        |
  1619.     addxl    d0,d0        |
  1620.     dbra    d4,1b        |
  1621. 2:
  1622. | Now round: we do it as follows: after the shifting we can write the
  1623. | fraction part as f + delta, where 1 < f < 2^25, and 0 <= delta <= 2.
  1624. | If delta < 1, do nothing. If delta > 1, add 1 to f. 
  1625. | If delta == 1, we make sure the rounded number will be even (odd?) 
  1626. | (after shifting).
  1627.     btst    IMM (0),d1    | is delta < 1?
  1628.     beq    2f        | if so, do not do anything
  1629.     orl    d2,d3        | is delta == 1?
  1630.     bne    1f        | if so round to even
  1631.     movel    d1,d3        | 
  1632.     andl    IMM (2),d3    | bit 1 is the last significant bit
  1633.     movel    IMM (0),d2    |
  1634.     addl    d3,d1        |
  1635.     addxl    d2,d0        |
  1636.     bra    2f        | 
  1637. 1:    movel    IMM (1),d3    | else add 1 
  1638.     movel    IMM (0),d2    |
  1639.     addl    d3,d1        |
  1640.     addxl    d2,d0
  1641. | Shift right once (because we used bit #DBL_MANT_DIG-32!).
  1642. 2:    lsrl    IMM (1),d0
  1643.     roxrl    IMM (1),d1        
  1644.  
  1645. | Now check again bit #DBL_MANT_DIG-32 (rounding could have produced a
  1646. | 'fraction overflow' ...).
  1647.     btst    IMM (DBL_MANT_DIG-32),d0    
  1648.     beq    1f
  1649.     lsrl    IMM (1),d0
  1650.     roxrl    IMM (1),d1
  1651.     addw    IMM (1),d4
  1652. 1:
  1653. | If bit #DBL_MANT_DIG-32-1 is clear we have a denormalized number, so we 
  1654. | have to put the exponent to zero and return a denormalized number.
  1655.     btst    IMM (DBL_MANT_DIG-32-1),d0
  1656.     beq    1f
  1657.     jmp    a0@
  1658. 1:    movel    IMM (0),d4
  1659.     jmp    a0@
  1660.  
  1661. Lround$to$zero:
  1662. Lround$to$plus:
  1663. Lround$to$minus:
  1664.     jmp    a0@
  1665. #endif /* L_double */
  1666.  
  1667. #ifdef  L_float
  1668.  
  1669.     .globl    SYM (_fpCCR)
  1670.     .globl  $_exception_handler
  1671.  
  1672. QUIET_NaN    = 0xffffffff
  1673. SIGNL_NaN    = 0x7f800001
  1674. INFINITY     = 0x7f800000
  1675.  
  1676. F_MAX_EXP      = 0xff
  1677. F_BIAS         = 126
  1678. FLT_MAX_EXP    = F_MAX_EXP - F_BIAS
  1679. FLT_MIN_EXP    = 1 - F_BIAS
  1680. FLT_MANT_DIG   = 24
  1681.  
  1682. INEXACT_RESULT         = 0x0001
  1683. UNDERFLOW         = 0x0002
  1684. OVERFLOW         = 0x0004
  1685. DIVIDE_BY_ZERO         = 0x0008
  1686. INVALID_OPERATION     = 0x0010
  1687.  
  1688. SINGLE_FLOAT = 1
  1689.  
  1690. NOOP         = 0
  1691. ADD          = 1
  1692. MULTIPLY     = 2
  1693. DIVIDE       = 3
  1694. NEGATE       = 4
  1695. COMPARE      = 5
  1696. EXTENDSFDF   = 6
  1697. TRUNCDFSF    = 7
  1698.  
  1699. UNKNOWN           = -1
  1700. ROUND_TO_NEAREST  = 0 | round result to nearest representable value
  1701. ROUND_TO_ZERO     = 1 | round result towards zero
  1702. ROUND_TO_PLUS     = 2 | round result towards plus infinity
  1703. ROUND_TO_MINUS    = 3 | round result towards minus infinity
  1704.  
  1705. | Entry points:
  1706.  
  1707.     .globl SYM (__addsf3)
  1708.     .globl SYM (__subsf3)
  1709.     .globl SYM (__mulsf3)
  1710.     .globl SYM (__divsf3)
  1711.     .globl SYM (__negsf2)
  1712.     .globl SYM (__cmpsf2)
  1713.  
  1714. | These are common routines to return and signal exceptions.    
  1715.  
  1716.     .text
  1717.     .even
  1718.  
  1719. Lf$den:
  1720. | Return and signal a denormalized number
  1721.     orl    d7,d0
  1722.     movew    IMM (UNDERFLOW),d7
  1723.     orw    IMM (INEXACT_RESULT),d7
  1724.     movew    IMM (SINGLE_FLOAT),d6
  1725.     jmp    $_exception_handler
  1726.  
  1727. Lf$infty:
  1728. Lf$overflow:
  1729. | Return a properly signed INFINITY and set the exception flags 
  1730.     movel    IMM (INFINITY),d0
  1731.     orl    d7,d0
  1732.     movew    IMM (OVERFLOW),d7
  1733.     orw    IMM (INEXACT_RESULT),d7
  1734.     movew    IMM (SINGLE_FLOAT),d6
  1735.     jmp    $_exception_handler
  1736.  
  1737. Lf$underflow:
  1738. | Return 0 and set the exception flags 
  1739.     movel    IMM (0),d0
  1740.     movew    IMM (UNDERFLOW),d7
  1741.     orw    IMM (INEXACT_RESULT),d7
  1742.     movew    IMM (SINGLE_FLOAT),d6
  1743.     jmp    $_exception_handler
  1744.  
  1745. Lf$inop:
  1746. | Return a quiet NaN and set the exception flags
  1747.     movel    IMM (QUIET_NaN),d0
  1748.     movew    IMM (INVALID_OPERATION),d7
  1749.     orw    IMM (INEXACT_RESULT),d7
  1750.     movew    IMM (SINGLE_FLOAT),d6
  1751.     jmp    $_exception_handler
  1752.  
  1753. Lf$div$0:
  1754. | Return a properly signed INFINITY and set the exception flags
  1755.     movel    IMM (INFINITY),d0
  1756.     orl    d7,d0
  1757.     movew    IMM (DIVIDE_BY_ZERO),d7
  1758.     orw    IMM (INEXACT_RESULT),d7
  1759.     movew    IMM (SINGLE_FLOAT),d6
  1760.     jmp    $_exception_handler
  1761.  
  1762. |=============================================================================
  1763. |=============================================================================
  1764. |                         single precision routines
  1765. |=============================================================================
  1766. |=============================================================================
  1767.  
  1768. | A single precision floating point number (float) has the format:
  1769. |
  1770. | struct _float {
  1771. |  unsigned int sign      : 1;  /* sign bit */ 
  1772. |  unsigned int exponent  : 8;  /* exponent, shifted by 126 */
  1773. |  unsigned int fraction  : 23; /* fraction */
  1774. | } float;
  1775. | 
  1776. | Thus sizeof(float) = 4 (32 bits). 
  1777. |
  1778. | All the routines are callable from C programs, and return the result 
  1779. | in the single register d0. They also preserve all registers except 
  1780. | d0-d1 and a0-a1.
  1781.  
  1782. |=============================================================================
  1783. |                              __subsf3
  1784. |=============================================================================
  1785.  
  1786. | float __subsf3(float, float);
  1787. SYM (__subsf3):
  1788.     bchg    IMM (31),sp@(8)    | change sign of second operand
  1789.                 | and fall through
  1790. |=============================================================================
  1791. |                              __addsf3
  1792. |=============================================================================
  1793.  
  1794. | float __addsf3(float, float);
  1795. SYM (__addsf3):
  1796.     link    a6,IMM (0)    | everything will be done in registers
  1797.     moveml    d2-d7,sp@-    | save all data registers but d0-d1
  1798.     movel    a6@(8),d0    | get first operand
  1799.     movel    a6@(12),d1    | get second operand
  1800.     movel    d0,d6        | get d0's sign bit '
  1801.     addl    d0,d0        | check and clear sign bit of a
  1802.     beq    Laddsf$b    | if zero return second operand
  1803.     movel    d1,d7        | save b's sign bit '
  1804.     addl    d1,d1        | get rid of sign bit
  1805.     beq    Laddsf$a    | if zero return first operand
  1806.  
  1807.     movel    d6,a0        | save signs in address registers
  1808.     movel    d7,a1        | so we can use d6 and d7
  1809.  
  1810. | Get the exponents and check for denormalized and/or infinity.
  1811.  
  1812.     movel    IMM (0x00ffffff),d4    | mask to get fraction
  1813.     movel    IMM (0x01000000),d5    | mask to put hidden bit back
  1814.  
  1815.     movel    d0,d6        | save a to get exponent
  1816.     andl    d4,d0        | get fraction in d0
  1817.     notl     d4        | make d4 into a mask for the exponent
  1818.     andl    d4,d6        | get exponent in d6
  1819.     beq    Laddsf$a$den    | branch if a is denormalized
  1820.     cmpl    d4,d6        | check for INFINITY or NaN
  1821.     beq    Laddsf$nf
  1822.     swap    d6        | put exponent into first word
  1823.     orl    d5,d0        | and put hidden bit back
  1824. Laddsf$1:
  1825. | Now we have a's exponent in d6 (second byte) and the mantissa in d0. '
  1826.     movel    d1,d7        | get exponent in d7
  1827.     andl    d4,d7        | 
  1828.     beq    Laddsf$b$den    | branch if b is denormalized
  1829.     cmpl    d4,d7        | check for INFINITY or NaN
  1830.     beq    Laddsf$nf
  1831.     swap    d7        | put exponent into first word
  1832.     notl     d4        | make d4 into a mask for the fraction
  1833.     andl    d4,d1        | get fraction in d1
  1834.     orl    d5,d1        | and put hidden bit back
  1835. Laddsf$2:
  1836. | Now we have b's exponent in d7 (second byte) and the mantissa in d1. '
  1837.  
  1838. | Note that the hidden bit corresponds to bit #FLT_MANT_DIG-1, and we 
  1839. | shifted right once, so bit #FLT_MANT_DIG is set (so we have one extra
  1840. | bit).
  1841.  
  1842.     movel    d1,d2        | move b to d2, since we want to use
  1843.                 | two registers to do the sum
  1844.     movel    IMM (0),d1    | and clear the new ones
  1845.     movel    d1,d3        |
  1846.  
  1847. | Here we shift the numbers in registers d0 and d1 so the exponents are the
  1848. | same, and put the largest exponent in d6. Note that we are using two
  1849. | registers for each number (see the discussion by D. Knuth in "Seminumerical 
  1850. | Algorithms").
  1851.     cmpw    d6,d7        | compare exponents
  1852.     beq    Laddsf$3    | if equal don't shift '
  1853.     bhi    5f        | branch if second exponent largest
  1854. 1:
  1855.     subl    d6,d7        | keep the largest exponent
  1856.     negl    d7
  1857.     lsrw    IMM (8),d7    | put difference in lower byte
  1858. | if difference is too large we don't shift (actually, we can just exit) '
  1859.     cmpw    IMM (FLT_MANT_DIG+2),d7        
  1860.     bge    Laddsf$b$small
  1861.     cmpw    IMM (16),d7    | if difference >= 16 swap
  1862.     bge    4f
  1863. 2:
  1864.     subw    IMM (1),d7
  1865. 3:    lsrl    IMM (1),d2    | shift right second operand
  1866.     roxrl    IMM (1),d3
  1867.     dbra    d7,3b
  1868.     bra    Laddsf$3
  1869. 4:
  1870.     movew    d2,d3
  1871.     swap    d3
  1872.     movew    d3,d2
  1873.     swap    d2
  1874.     subw    IMM (16),d7
  1875.     bne    2b        | if still more bits, go back to normal case
  1876.     bra    Laddsf$3
  1877. 5:
  1878.     exg    d6,d7        | exchange the exponents
  1879.     subl    d6,d7        | keep the largest exponent
  1880.     negl    d7        |
  1881.     lsrw    IMM (8),d7    | put difference in lower byte
  1882. | if difference is too large we don't shift (and exit!) '
  1883.     cmpw    IMM (FLT_MANT_DIG+2),d7        
  1884.     bge    Laddsf$a$small
  1885.     cmpw    IMM (16),d7    | if difference >= 16 swap
  1886.     bge    8f
  1887. 6:
  1888.     subw    IMM (1),d7
  1889. 7:    lsrl    IMM (1),d0    | shift right first operand
  1890.     roxrl    IMM (1),d1
  1891.     dbra    d7,7b
  1892.     bra    Laddsf$3
  1893. 8:
  1894.     movew    d0,d1
  1895.     swap    d1
  1896.     movew    d1,d0
  1897.     swap    d0
  1898.     subw    IMM (16),d7
  1899.     bne    6b        | if still more bits, go back to normal case
  1900.                 | otherwise we fall through
  1901.  
  1902. | Now we have a in d0-d1, b in d2-d3, and the largest exponent in d6 (the
  1903. | signs are stored in a0 and a1).
  1904.  
  1905. Laddsf$3:
  1906. | Here we have to decide whether to add or subtract the numbers
  1907.     exg    d6,a0        | get signs back
  1908.     exg    d7,a1        | and save the exponents
  1909.     eorl    d6,d7        | combine sign bits
  1910.     bmi    Lsubsf$0    | if negative a and b have opposite 
  1911.                 | sign so we actually subtract the
  1912.                 | numbers
  1913.  
  1914. | Here we have both positive or both negative
  1915.     exg    d6,a0        | now we have the exponent in d6
  1916.     movel    a0,d7        | and sign in d7
  1917.     andl    IMM (0x80000000),d7
  1918. | Here we do the addition.
  1919.     addl    d3,d1
  1920.     addxl    d2,d0
  1921. | Note: now we have d2, d3, d4 and d5 to play with! 
  1922.  
  1923. | Put the exponent, in the first byte, in d2, to use the "standard" rounding
  1924. | routines:
  1925.     movel    d6,d2
  1926.     lsrw    IMM (8),d2
  1927.  
  1928. | Before rounding normalize so bit #FLT_MANT_DIG is set (we will consider
  1929. | the case of denormalized numbers in the rounding routine itself).
  1930. | As in the addition (not in the subtraction!) we could have set 
  1931. | one more bit we check this:
  1932.     btst    IMM (FLT_MANT_DIG+1),d0    
  1933.     beq    1f
  1934.     lsrl    IMM (1),d0
  1935.     roxrl    IMM (1),d1
  1936.     addl    IMM (1),d2
  1937. 1:
  1938.     lea    Laddsf$4,a0    | to return from rounding routine
  1939.     lea    SYM (_fpCCR),a1    | check the rounding mode
  1940.     movew    a1@(6),d6    | rounding mode in d6
  1941.     beq    Lround$to$nearest
  1942.     cmpw    IMM (ROUND_TO_PLUS),d6
  1943.     bhi    Lround$to$minus
  1944.     blt    Lround$to$zero
  1945.     bra    Lround$to$plus
  1946. Laddsf$4:
  1947. | Put back the exponent, but check for overflow.
  1948.     cmpw    IMM (0xff),d2
  1949.     bhi    1f
  1950.     bclr    IMM (FLT_MANT_DIG-1),d0
  1951.     lslw    IMM (7),d2
  1952.     swap    d2
  1953.     orl    d2,d0
  1954.     bra    Laddsf$ret
  1955. 1:
  1956.     movew    IMM (ADD),d5
  1957.     bra    Lf$overflow
  1958.  
  1959. Lsubsf$0:
  1960. | We are here if a > 0 and b < 0 (sign bits cleared).
  1961. | Here we do the subtraction.
  1962.     movel    d6,d7        | put sign in d7
  1963.     andl    IMM (0x80000000),d7
  1964.  
  1965.     subl    d3,d1        | result in d0-d1
  1966.     subxl    d2,d0        |
  1967.     beq    Laddsf$ret    | if zero just exit
  1968.     bpl    1f        | if positive skip the following
  1969.     bchg    IMM (31),d7    | change sign bit in d7
  1970.     negl    d1
  1971.     negxl    d0
  1972. 1:
  1973.     exg    d2,a0        | now we have the exponent in d2
  1974.     lsrw    IMM (8),d2    | put it in the first byte
  1975.  
  1976. | Now d0-d1 is positive and the sign bit is in d7.
  1977.  
  1978. | Note that we do not have to normalize, since in the subtraction bit
  1979. | #FLT_MANT_DIG+1 is never set, and denormalized numbers are handled by
  1980. | the rounding routines themselves.
  1981.     lea    Lsubsf$1,a0    | to return from rounding routine
  1982.     lea    SYM (_fpCCR),a1    | check the rounding mode
  1983.     movew    a1@(6),d6    | rounding mode in d6
  1984.     beq    Lround$to$nearest
  1985.     cmpw    IMM (ROUND_TO_PLUS),d6
  1986.     bhi    Lround$to$minus
  1987.     blt    Lround$to$zero
  1988.     bra    Lround$to$plus
  1989. Lsubsf$1:
  1990. | Put back the exponent (we can't have overflow!). '
  1991.     bclr    IMM (FLT_MANT_DIG-1),d0
  1992.     lslw    IMM (7),d2
  1993.     swap    d2
  1994.     orl    d2,d0
  1995.     bra    Laddsf$ret
  1996.  
  1997. | If one of the numbers was too small (difference of exponents >= 
  1998. | FLT_MANT_DIG+2) we return the other (and now we don't have to '
  1999. | check for finiteness or zero).
  2000. Laddsf$a$small:
  2001.     movel    a6@(12),d0
  2002.     lea    SYM (_fpCCR),a0
  2003.     movew    IMM (0),a0@
  2004.     moveml    sp@+,d2-d7    | restore data registers
  2005.     unlk    a6        | and return
  2006.     rts
  2007.  
  2008. Laddsf$b$small:
  2009.     movel    a6@(8),d0
  2010.     lea    SYM (_fpCCR),a0
  2011.     movew    IMM (0),a0@
  2012.     moveml    sp@+,d2-d7    | restore data registers
  2013.     unlk    a6        | and return
  2014.     rts
  2015.  
  2016. | If the numbers are denormalized remember to put exponent equal to 1.
  2017.  
  2018. Laddsf$a$den:
  2019.     movel    d5,d6        | d5 contains 0x01000000
  2020.     swap    d6
  2021.     bra    Laddsf$1
  2022.  
  2023. Laddsf$b$den:
  2024.     movel    d5,d7
  2025.     swap    d7
  2026.     notl     d4        | make d4 into a mask for the fraction
  2027.                 | (this was not executed after the jump)
  2028.     bra    Laddsf$2
  2029.  
  2030. | The rest is mainly code for the different results which can be 
  2031. | returned (checking always for +/-INFINITY and NaN).
  2032.  
  2033. Laddsf$b:
  2034. | Return b (if a is zero).
  2035.     movel    a6@(12),d0
  2036.     bra    1f
  2037. Laddsf$a:
  2038. | Return a (if b is zero).
  2039.     movel    a6@(8),d0
  2040. 1:
  2041.     movew    IMM (ADD),d5
  2042. | We have to check for NaN and +/-infty.
  2043.     movel    d0,d7
  2044.     andl    IMM (0x80000000),d7    | put sign in d7
  2045.     bclr    IMM (31),d0        | clear sign
  2046.     cmpl    IMM (INFINITY),d0    | check for infty or NaN
  2047.     bge    2f
  2048.     movel    d0,d0        | check for zero (we do this because we don't '
  2049.     bne    Laddsf$ret    | want to return -0 by mistake
  2050.     bclr    IMM (31),d7    | if zero be sure to clear sign
  2051.     bra    Laddsf$ret    | if everything OK just return
  2052. 2:
  2053. | The value to be returned is either +/-infty or NaN
  2054.     andl    IMM (0x007fffff),d0    | check for NaN
  2055.     bne    Lf$inop            | if mantissa not zero is NaN
  2056.     bra    Lf$infty
  2057.  
  2058. Laddsf$ret:
  2059. | Normal exit (a and b nonzero, result is not NaN nor +/-infty).
  2060. | We have to clear the exception flags (just the exception type).
  2061.     lea    SYM (_fpCCR),a0
  2062.     movew    IMM (0),a0@
  2063.     orl    d7,d0        | put sign bit
  2064.     moveml    sp@+,d2-d7    | restore data registers
  2065.     unlk    a6        | and return
  2066.     rts
  2067.  
  2068. Laddsf$ret$den:
  2069. | Return a denormalized number (for addition we don't signal underflow) '
  2070.     lsrl    IMM (1),d0    | remember to shift right back once
  2071.     bra    Laddsf$ret    | and return
  2072.  
  2073. | Note: when adding two floats of the same sign if either one is 
  2074. | NaN we return NaN without regard to whether the other is finite or 
  2075. | not. When subtracting them (i.e., when adding two numbers of 
  2076. | opposite signs) things are more complicated: if both are INFINITY 
  2077. | we return NaN, if only one is INFINITY and the other is NaN we return
  2078. | NaN, but if it is finite we return INFINITY with the corresponding sign.
  2079.  
  2080. Laddsf$nf:
  2081.     movew    IMM (ADD),d5
  2082. | This could be faster but it is not worth the effort, since it is not
  2083. | executed very often. We sacrifice speed for clarity here.
  2084.     movel    a6@(8),d0    | get the numbers back (remember that we
  2085.     movel    a6@(12),d1    | did some processing already)
  2086.     movel    IMM (INFINITY),d4 | useful constant (INFINITY)
  2087.     movel    d0,d2        | save sign bits
  2088.     movel    d1,d3
  2089.     bclr    IMM (31),d0    | clear sign bits
  2090.     bclr    IMM (31),d1
  2091. | We know that one of them is either NaN of +/-INFINITY
  2092. | Check for NaN (if either one is NaN return NaN)
  2093.     cmpl    d4,d0        | check first a (d0)
  2094.     bhi    Lf$inop        
  2095.     cmpl    d4,d1        | check now b (d1)
  2096.     bhi    Lf$inop        
  2097. | Now comes the check for +/-INFINITY. We know that both are (maybe not
  2098. | finite) numbers, but we have to check if both are infinite whether we
  2099. | are adding or subtracting them.
  2100.     eorl    d3,d2        | to check sign bits
  2101.     bmi    1f
  2102.     movel    d0,d7
  2103.     andl    IMM (0x80000000),d7    | get (common) sign bit
  2104.     bra    Lf$infty
  2105. 1:
  2106. | We know one (or both) are infinite, so we test for equality between the
  2107. | two numbers (if they are equal they have to be infinite both, so we
  2108. | return NaN).
  2109.     cmpl    d1,d0        | are both infinite?
  2110.     beq    Lf$inop        | if so return NaN
  2111.  
  2112.     movel    d0,d7
  2113.     andl    IMM (0x80000000),d7 | get a's sign bit '
  2114.     cmpl    d4,d0        | test now for infinity
  2115.     beq    Lf$infty    | if a is INFINITY return with this sign
  2116.     bchg    IMM (31),d7    | else we know b is INFINITY and has
  2117.     bra    Lf$infty    | the opposite sign
  2118.  
  2119. |=============================================================================
  2120. |                             __mulsf3
  2121. |=============================================================================
  2122.  
  2123. | float __mulsf3(float, float);
  2124. SYM (__mulsf3):
  2125.     link    a6,IMM (0)
  2126.     moveml    d2-d7,sp@-
  2127.     movel    a6@(8),d0    | get a into d0
  2128.     movel    a6@(12),d1    | and b into d1
  2129.     movel    d0,d7        | d7 will hold the sign of the product
  2130.     eorl    d1,d7        |
  2131.     andl    IMM (0x80000000),d7
  2132.     movel    IMM (INFINITY),d6    | useful constant (+INFINITY)
  2133.     movel    d6,d5            | another (mask for fraction)
  2134.     notl    d5            |
  2135.     movel    IMM (0x00800000),d4    | this is to put hidden bit back
  2136.     bclr    IMM (31),d0        | get rid of a's sign bit '
  2137.     movel    d0,d2            |
  2138.     beq    Lmulsf$a$0        | branch if a is zero
  2139.     bclr    IMM (31),d1        | get rid of b's sign bit '
  2140.     movel    d1,d3        |
  2141.     beq    Lmulsf$b$0    | branch if b is zero
  2142.     cmpl    d6,d0        | is a big?
  2143.     bhi    Lmulsf$inop    | if a is NaN return NaN
  2144.     beq    Lmulsf$inf    | if a is INFINITY we have to check b
  2145.     cmpl    d6,d1        | now compare b with INFINITY
  2146.     bhi    Lmulsf$inop    | is b NaN?
  2147.     beq    Lmulsf$overflow | is b INFINITY?
  2148. | Here we have both numbers finite and nonzero (and with no sign bit).
  2149. | Now we get the exponents into d2 and d3.
  2150.     andl    d6,d2        | and isolate exponent in d2
  2151.     beq    Lmulsf$a$den    | if exponent is zero we have a denormalized
  2152.     andl    d5,d0        | and isolate fraction
  2153.     orl    d4,d0        | and put hidden bit back
  2154.     swap    d2        | I like exponents in the first byte
  2155.     lsrw    IMM (7),d2    | 
  2156. Lmulsf$1:            | number
  2157.     andl    d6,d3        |
  2158.     beq    Lmulsf$b$den    |
  2159.     andl    d5,d1        |
  2160.     orl    d4,d1        |
  2161.     swap    d3        |
  2162.     lsrw    IMM (7),d3    |
  2163. Lmulsf$2:            |
  2164.     addw    d3,d2        | add exponents
  2165.     subw    IMM (F_BIAS+1),d2 | and subtract bias (plus one)
  2166.  
  2167. | We are now ready to do the multiplication. The situation is as follows:
  2168. | both a and b have bit FLT_MANT_DIG-1 set (even if they were 
  2169. | denormalized to start with!), which means that in the product 
  2170. | bit 2*(FLT_MANT_DIG-1) (that is, bit 2*FLT_MANT_DIG-2-32 of the 
  2171. | high long) is set. 
  2172.  
  2173. | To do the multiplication let us move the number a little bit around ...
  2174.     movel    d1,d6        | second operand in d6
  2175.     movel    d0,d5        | first operand in d4-d5
  2176.     movel    IMM (0),d4
  2177.     movel    d4,d1        | the sums will go in d0-d1
  2178.     movel    d4,d0
  2179.  
  2180. | now bit FLT_MANT_DIG-1 becomes bit 31:
  2181.     lsll    IMM (31-FLT_MANT_DIG+1),d6        
  2182.  
  2183. | Start the loop (we loop #FLT_MANT_DIG times):
  2184.     movew    IMM (FLT_MANT_DIG-1),d3    
  2185. 1:    addl    d1,d1        | shift sum 
  2186.     addxl    d0,d0
  2187.     lsll    IMM (1),d6    | get bit bn
  2188.     bcc    2f        | if not set skip sum
  2189.     addl    d5,d1        | add a
  2190.     addxl    d4,d0
  2191. 2:    dbf    d3,1b        | loop back
  2192.  
  2193. | Now we have the product in d0-d1, with bit (FLT_MANT_DIG - 1) + FLT_MANT_DIG
  2194. | (mod 32) of d0 set. The first thing to do now is to normalize it so bit 
  2195. | FLT_MANT_DIG is set (to do the rounding).
  2196.     rorl    IMM (6),d1
  2197.     swap    d1
  2198.     movew    d1,d3
  2199.     andw    IMM (0x03ff),d3
  2200.     andw    IMM (0xfd00),d1
  2201.     lsll    IMM (8),d0
  2202.     addl    d0,d0
  2203.     addl    d0,d0
  2204.     orw    d3,d0
  2205.  
  2206.     movew    IMM (MULTIPLY),d5
  2207.     
  2208.     btst    IMM (FLT_MANT_DIG+1),d0
  2209.     beq    Lround$exit
  2210.     lsrl    IMM (1),d0
  2211.     roxrl    IMM (1),d1
  2212.     addw    IMM (1),d2
  2213.     bra    Lround$exit
  2214.  
  2215. Lmulsf$inop:
  2216.     movew    IMM (MULTIPLY),d5
  2217.     bra    Lf$inop
  2218.  
  2219. Lmulsf$overflow:
  2220.     movew    IMM (MULTIPLY),d5
  2221.     bra    Lf$overflow
  2222.  
  2223. Lmulsf$inf:
  2224.     movew    IMM (MULTIPLY),d5
  2225. | If either is NaN return NaN; else both are (maybe infinite) numbers, so
  2226. | return INFINITY with the correct sign (which is in d7).
  2227.     cmpl    d6,d1        | is b NaN?
  2228.     bhi    Lf$inop        | if so return NaN
  2229.     bra    Lf$overflow    | else return +/-INFINITY
  2230.  
  2231. | If either number is zero return zero, unless the other is +/-INFINITY, 
  2232. | or NaN, in which case we return NaN.
  2233. Lmulsf$b$0:
  2234. | Here d1 (==b) is zero.
  2235.     movel    d1,d0        | put b into d0 (just a zero)
  2236.     movel    a6@(8),d1    | get a again to check for non-finiteness
  2237.     bra    1f
  2238. Lmulsf$a$0:
  2239.     movel    a6@(12),d1    | get b again to check for non-finiteness
  2240. 1:    bclr    IMM (31),d1    | clear sign bit 
  2241.     cmpl    IMM (INFINITY),d1 | and check for a large exponent
  2242.     bge    Lf$inop        | if b is +/-INFINITY or NaN return NaN
  2243.     lea    SYM (_fpCCR),a0    | else return zero
  2244.     movew    IMM (0),a0@    | 
  2245.     moveml    sp@+,d2-d7    | 
  2246.     unlk    a6        | 
  2247.     rts            | 
  2248.  
  2249. | If a number is denormalized we put an exponent of 1 but do not put the 
  2250. | hidden bit back into the fraction; instead we shift left until bit 23
  2251. | (the hidden bit) is set, adjusting the exponent accordingly. We do this
  2252. | to ensure that the product of the fractions is close to 1.
  2253. Lmulsf$a$den:
  2254.     movel    IMM (1),d2
  2255.     andl    d5,d0
  2256. 1:    addl    d0,d0        | shift a left (until bit 23 is set)
  2257.     subw    IMM (1),d2    | and adjust exponent
  2258.     btst    IMM (FLT_MANT_DIG-1),d0
  2259.     bne    Lmulsf$1    |
  2260.     bra    1b        | else loop back
  2261.  
  2262. Lmulsf$b$den:
  2263.     movel    IMM (1),d3
  2264.     andl    d5,d1
  2265. 1:    addl    d1,d1        | shift b left until bit 23 is set
  2266.     subw    IMM (1),d3    | and adjust exponent
  2267.     btst    IMM (FLT_MANT_DIG-1),d1
  2268.     bne    Lmulsf$2    |
  2269.     bra    1b        | else loop back
  2270.  
  2271. |=============================================================================
  2272. |                             __divsf3
  2273. |=============================================================================
  2274.  
  2275. | float __divsf3(float, float);
  2276. SYM (__divsf3):
  2277.     link    a6,IMM (0)
  2278.     moveml    d2-d7,sp@-
  2279.     movel    a6@(8),d0        | get a into d0
  2280.     movel    a6@(12),d1        | and b into d1
  2281.     movel    d0,d7            | d7 will hold the sign of the result
  2282.     eorl    d1,d7            |
  2283.     andl    IMM (0x80000000),d7    | 
  2284.     movel    IMM (INFINITY),d6    | useful constant (+INFINITY)
  2285.     movel    d6,d5            | another (mask for fraction)
  2286.     notl    d5            |
  2287.     movel    IMM (0x00800000),d4    | this is to put hidden bit back
  2288.     bclr    IMM (31),d0        | get rid of a's sign bit '
  2289.     movel    d0,d2            |
  2290.     beq    Ldivsf$a$0        | branch if a is zero
  2291.     bclr    IMM (31),d1        | get rid of b's sign bit '
  2292.     movel    d1,d3            |
  2293.     beq    Ldivsf$b$0        | branch if b is zero
  2294.     cmpl    d6,d0            | is a big?
  2295.     bhi    Ldivsf$inop        | if a is NaN return NaN
  2296.     beq    Ldivsf$inf        | if a is INFINITY we have to check b
  2297.     cmpl    d6,d1            | now compare b with INFINITY 
  2298.     bhi    Ldivsf$inop        | if b is NaN return NaN
  2299.     beq    Ldivsf$underflow
  2300. | Here we have both numbers finite and nonzero (and with no sign bit).
  2301. | Now we get the exponents into d2 and d3 and normalize the numbers to
  2302. | ensure that the ratio of the fractions is close to 1. We do this by
  2303. | making sure that bit #FLT_MANT_DIG-1 (hidden bit) is set.
  2304.     andl    d6,d2        | and isolate exponent in d2
  2305.     beq    Ldivsf$a$den    | if exponent is zero we have a denormalized
  2306.     andl    d5,d0        | and isolate fraction
  2307.     orl    d4,d0        | and put hidden bit back
  2308.     swap    d2        | I like exponents in the first byte
  2309.     lsrw    IMM (7),d2    | 
  2310. Ldivsf$1:            | 
  2311.     andl    d6,d3        |
  2312.     beq    Ldivsf$b$den    |
  2313.     andl    d5,d1        |
  2314.     orl    d4,d1        |
  2315.     swap    d3        |
  2316.     lsrw    IMM (7),d3    |
  2317. Ldivsf$2:            |
  2318.     subw    d3,d2        | subtract exponents
  2319.      addw    IMM (F_BIAS),d2    | and add bias
  2320.  
  2321. | We are now ready to do the division. We have prepared things in such a way
  2322. | that the ratio of the fractions will be less than 2 but greater than 1/2.
  2323. | At this point the registers in use are:
  2324. | d0    holds a (first operand, bit FLT_MANT_DIG=0, bit FLT_MANT_DIG-1=1)
  2325. | d1    holds b (second operand, bit FLT_MANT_DIG=1)
  2326. | d2    holds the difference of the exponents, corrected by the bias
  2327. | d7    holds the sign of the ratio
  2328. | d4, d5, d6 hold some constants
  2329.     movel    d7,a0        | d6-d7 will hold the ratio of the fractions
  2330.     movel    IMM (0),d6    | 
  2331.     movel    d6,d7
  2332.  
  2333.     movew    IMM (FLT_MANT_DIG+1),d3
  2334. 1:    cmpl    d0,d1        | is a < b?
  2335.     bhi    2f        |
  2336.     bset    d3,d6        | set a bit in d6
  2337.     subl    d1,d0        | if a >= b  a <-- a-b
  2338.     beq    3f        | if a is zero, exit
  2339. 2:    addl    d0,d0        | multiply a by 2
  2340.     dbra    d3,1b
  2341.  
  2342. | Now we keep going to set the sticky bit ...
  2343.     movew    IMM (FLT_MANT_DIG),d3
  2344. 1:    cmpl    d0,d1
  2345.     ble    2f
  2346.     addl    d0,d0
  2347.     dbra    d3,1b
  2348.     movel    IMM (0),d1
  2349.     bra    3f
  2350. 2:    movel    IMM (0),d1
  2351.     subw    IMM (FLT_MANT_DIG),d3
  2352.     addw    IMM (31),d3
  2353.     bset    d3,d1
  2354. 3:
  2355.     movel    d6,d0        | put the ratio in d0-d1
  2356.     movel    a0,d7        | get sign back
  2357.  
  2358. | Because of the normalization we did before we are guaranteed that 
  2359. | d0 is smaller than 2^26 but larger than 2^24. Thus bit 26 is not set,
  2360. | bit 25 could be set, and if it is not set then bit 24 is necessarily set.
  2361.     btst    IMM (FLT_MANT_DIG+1),d0        
  2362.     beq    1f              | if it is not set, then bit 24 is set
  2363.     lsrl    IMM (1),d0    |
  2364.     addw    IMM (1),d2    |
  2365. 1:
  2366. | Now round, check for over- and underflow, and exit.
  2367.     movew    IMM (DIVIDE),d5
  2368.     bra    Lround$exit
  2369.  
  2370. Ldivsf$inop:
  2371.     movew    IMM (DIVIDE),d5
  2372.     bra    Lf$inop
  2373.  
  2374. Ldivsf$overflow:
  2375.     movew    IMM (DIVIDE),d5
  2376.     bra    Lf$overflow
  2377.  
  2378. Ldivsf$underflow:
  2379.     movew    IMM (DIVIDE),d5
  2380.     bra    Lf$underflow
  2381.  
  2382. Ldivsf$a$0:
  2383.     movew    IMM (DIVIDE),d5
  2384. | If a is zero check to see whether b is zero also. In that case return
  2385. | NaN; then check if b is NaN, and return NaN also in that case. Else
  2386. | return zero.
  2387.     andl    IMM (0x7fffffff),d1    | clear sign bit and test b
  2388.     beq    Lf$inop            | if b is also zero return NaN
  2389.     cmpl    IMM (INFINITY),d1    | check for NaN
  2390.     bhi    Lf$inop            | 
  2391.     movel    IMM (0),d0        | else return zero
  2392.     lea    SYM (_fpCCR),a0        |
  2393.     movew    IMM (0),a0@        |
  2394.     moveml    sp@+,d2-d7        | 
  2395.     unlk    a6            | 
  2396.     rts                | 
  2397.     
  2398. Ldivsf$b$0:
  2399.     movew    IMM (DIVIDE),d5
  2400. | If we got here a is not zero. Check if a is NaN; in that case return NaN,
  2401. | else return +/-INFINITY. Remember that a is in d0 with the sign bit 
  2402. | cleared already.
  2403.     cmpl    IMM (INFINITY),d0    | compare d0 with INFINITY
  2404.     bhi    Lf$inop            | if larger it is NaN
  2405.     bra    Lf$div$0        | else signal DIVIDE_BY_ZERO
  2406.  
  2407. Ldivsf$inf:
  2408.     movew    IMM (DIVIDE),d5
  2409. | If a is INFINITY we have to check b
  2410.     cmpl    IMM (INFINITY),d1    | compare b with INFINITY 
  2411.     bge    Lf$inop            | if b is NaN or INFINITY return NaN
  2412.     bra    Lf$overflow        | else return overflow
  2413.  
  2414. | If a number is denormalized we put an exponent of 1 but do not put the 
  2415. | bit back into the fraction.
  2416. Ldivsf$a$den:
  2417.     movel    IMM (1),d2
  2418.     andl    d5,d0
  2419. 1:    addl    d0,d0        | shift a left until bit FLT_MANT_DIG-1 is set
  2420.     subw    IMM (1),d2    | and adjust exponent
  2421.     btst    IMM (FLT_MANT_DIG-1),d0
  2422.     bne    Ldivsf$1
  2423.     bra    1b
  2424.  
  2425. Ldivsf$b$den:
  2426.     movel    IMM (1),d3
  2427.     andl    d5,d1
  2428. 1:    addl    d1,d1        | shift b left until bit FLT_MANT_DIG is set
  2429.     subw    IMM (1),d3    | and adjust exponent
  2430.     btst    IMM (FLT_MANT_DIG-1),d1
  2431.     bne    Ldivsf$2
  2432.     bra    1b
  2433.  
  2434. Lround$exit:
  2435. | This is a common exit point for __mulsf3 and __divsf3. 
  2436.  
  2437. | First check for underlow in the exponent:
  2438.     cmpw    IMM (-FLT_MANT_DIG-1),d2        
  2439.     blt    Lf$underflow    
  2440. | It could happen that the exponent is less than 1, in which case the 
  2441. | number is denormalized. In this case we shift right and adjust the 
  2442. | exponent until it becomes 1 or the fraction is zero (in the latter case 
  2443. | we signal underflow and return zero).
  2444.     movel    IMM (0),d6    | d6 is used temporarily
  2445.     cmpw    IMM (1),d2    | if the exponent is less than 1 we 
  2446.     bge    2f        | have to shift right (denormalize)
  2447. 1:    addw    IMM (1),d2    | adjust the exponent
  2448.     lsrl    IMM (1),d0    | shift right once 
  2449.     roxrl    IMM (1),d1    |
  2450.     roxrl    IMM (1),d6    | d6 collect bits we would lose otherwise
  2451.     cmpw    IMM (1),d2    | is the exponent 1 already?
  2452.     beq    2f        | if not loop back
  2453.     bra    1b              |
  2454.     bra    Lf$underflow    | safety check, shouldn't execute '
  2455. 2:    orl    d6,d1        | this is a trick so we don't lose  '
  2456.                 | the extra bits which were flushed right
  2457. | Now call the rounding routine (which takes care of denormalized numbers):
  2458.     lea    Lround$0,a0    | to return from rounding routine
  2459.     lea    SYM (_fpCCR),a1    | check the rounding mode
  2460.     movew    a1@(6),d6    | rounding mode in d6
  2461.     beq    Lround$to$nearest
  2462.     cmpw    IMM (ROUND_TO_PLUS),d6
  2463.     bhi    Lround$to$minus
  2464.     blt    Lround$to$zero
  2465.     bra    Lround$to$plus
  2466. Lround$0:
  2467. | Here we have a correctly rounded result (either normalized or denormalized).
  2468.  
  2469. | Here we should have either a normalized number or a denormalized one, and
  2470. | the exponent is necessarily larger or equal to 1 (so we don't have to  '
  2471. | check again for underflow!). We have to check for overflow or for a 
  2472. | denormalized number (which also signals underflow).
  2473. | Check for overflow (i.e., exponent >= 255).
  2474.     cmpw    IMM (0x00ff),d2
  2475.     bge    Lf$overflow
  2476. | Now check for a denormalized number (exponent==0).
  2477.     movew    d2,d2
  2478.     beq    Lf$den
  2479. 1:
  2480. | Put back the exponents and sign and return.
  2481.     lslw    IMM (7),d2    | exponent back to fourth byte
  2482.     bclr    IMM (FLT_MANT_DIG-1),d0
  2483.     swap    d0        | and put back exponent
  2484.     orw    d2,d0        | 
  2485.     swap    d0        |
  2486.     orl    d7,d0        | and sign also
  2487.  
  2488.     lea    SYM (_fpCCR),a0
  2489.     movew    IMM (0),a0@
  2490.     moveml    sp@+,d2-d7
  2491.     unlk    a6
  2492.     rts
  2493.  
  2494. |=============================================================================
  2495. |                             __negsf2
  2496. |=============================================================================
  2497.  
  2498. | This is trivial and could be shorter if we didn't bother checking for NaN '
  2499. | and +/-INFINITY.
  2500.  
  2501. | float __negsf2(float);
  2502. SYM (__negsf2):
  2503.     link    a6,IMM (0)
  2504.     moveml    d2-d7,sp@-
  2505.     movew    IMM (NEGATE),d5
  2506.     movel    a6@(8),d0    | get number to negate in d0
  2507.     bchg    IMM (31),d0    | negate
  2508.     movel    d0,d1        | make a positive copy
  2509.     bclr    IMM (31),d1    |
  2510.     tstl    d1        | check for zero
  2511.     beq    2f        | if zero (either sign) return +zero
  2512.     cmpl    IMM (INFINITY),d1 | compare to +INFINITY
  2513.     blt    1f        |
  2514.     bhi    Lf$inop        | if larger (fraction not zero) is NaN
  2515.     movel    d0,d7        | else get sign and return INFINITY
  2516.     andl    IMM (0x80000000),d7
  2517.     bra    Lf$infty        
  2518. 1:    lea    SYM (_fpCCR),a0
  2519.     movew    IMM (0),a0@
  2520.     moveml    sp@+,d2-d7
  2521.     unlk    a6
  2522.     rts
  2523. 2:    bclr    IMM (31),d0
  2524.     bra    1b
  2525.  
  2526. |=============================================================================
  2527. |                             __cmpsf2
  2528. |=============================================================================
  2529.  
  2530. GREATER =  1
  2531. LESS    = -1
  2532. EQUAL   =  0
  2533.  
  2534. | int __cmpsf2(float, float);
  2535. SYM (__cmpsf2):
  2536.     link    a6,IMM (0)
  2537.     moveml    d2-d7,sp@-     | save registers
  2538.     movew    IMM (COMPARE),d5
  2539.     movel    a6@(8),d0    | get first operand
  2540.     movel    a6@(12),d1    | get second operand
  2541. | Check if either is NaN, and in that case return garbage and signal
  2542. | INVALID_OPERATION. Check also if either is zero, and clear the signs
  2543. | if necessary.
  2544.     movel    d0,d6
  2545.     andl    IMM (0x7fffffff),d0
  2546.     beq    Lcmpsf$a$0
  2547.     cmpl    IMM (0x7f800000),d0
  2548.     bhi    Lf$inop
  2549. Lcmpsf$1:
  2550.     movel    d1,d7
  2551.     andl    IMM (0x7fffffff),d1
  2552.     beq    Lcmpsf$b$0
  2553.     cmpl    IMM (0x7f800000),d1
  2554.     bhi    Lf$inop
  2555. Lcmpsf$2:
  2556. | Check the signs
  2557.     eorl    d6,d7
  2558.     bpl    1f
  2559. | If the signs are not equal check if a >= 0
  2560.     tstl    d6
  2561.     bpl    Lcmpsf$a$gt$b    | if (a >= 0 && b < 0) => a > b
  2562.     bmi    Lcmpsf$b$gt$a    | if (a < 0 && b >= 0) => a < b
  2563. 1:
  2564. | If the signs are equal check for < 0
  2565.     tstl    d6
  2566.     bpl    1f
  2567. | If both are negative exchange them
  2568.     exg    d0,d1
  2569. 1:
  2570. | Now that they are positive we just compare them as longs (does this also
  2571. | work for denormalized numbers?).
  2572.     cmpl    d0,d1
  2573.     bhi    Lcmpsf$b$gt$a    | |b| > |a|
  2574.     bne    Lcmpsf$a$gt$b    | |b| < |a|
  2575. | If we got here a == b.
  2576.     movel    IMM (EQUAL),d0
  2577.     moveml    sp@+,d2-d7     | put back the registers
  2578.     unlk    a6
  2579.     rts
  2580. Lcmpsf$a$gt$b:
  2581.     movel    IMM (GREATER),d0
  2582.     moveml    sp@+,d2-d7     | put back the registers
  2583.     unlk    a6
  2584.     rts
  2585. Lcmpsf$b$gt$a:
  2586.     movel    IMM (LESS),d0
  2587.     moveml    sp@+,d2-d7     | put back the registers
  2588.     unlk    a6
  2589.     rts
  2590.  
  2591. Lcmpsf$a$0:    
  2592.     bclr    IMM (31),d6
  2593.     bra    Lcmpsf$1
  2594. Lcmpsf$b$0:
  2595.     bclr    IMM (31),d7
  2596.     bra    Lcmpsf$2
  2597.  
  2598. |=============================================================================
  2599. |                           rounding routines
  2600. |=============================================================================
  2601.  
  2602. | The rounding routines expect the number to be normalized in registers
  2603. | d0-d1, with the exponent in register d2. They assume that the 
  2604. | exponent is larger or equal to 1. They return a properly normalized number
  2605. | if possible, and a denormalized number otherwise. The exponent is returned
  2606. | in d2.
  2607.  
  2608. Lround$to$nearest:
  2609. | We now normalize as suggested by D. Knuth ("Seminumerical Algorithms"):
  2610. | Here we assume that the exponent is not too small (this should be checked
  2611. | before entering the rounding routine), but the number could be denormalized.
  2612.  
  2613. | Check for denormalized numbers:
  2614. 1:    btst    IMM (FLT_MANT_DIG),d0
  2615.     bne    2f        | if set the number is normalized
  2616. | Normalize shifting left until bit #FLT_MANT_DIG is set or the exponent 
  2617. | is one (remember that a denormalized number corresponds to an 
  2618. | exponent of -F_BIAS+1).
  2619.     cmpw    IMM (1),d2    | remember that the exponent is at least one
  2620.      beq    2f        | an exponent of one means denormalized
  2621.     addl    d1,d1        | else shift and adjust the exponent
  2622.     addxl    d0,d0        |
  2623.     dbra    d2,1b        |
  2624. 2:
  2625. | Now round: we do it as follows: after the shifting we can write the
  2626. | fraction part as f + delta, where 1 < f < 2^25, and 0 <= delta <= 2.
  2627. | If delta < 1, do nothing. If delta > 1, add 1 to f. 
  2628. | If delta == 1, we make sure the rounded number will be even (odd?) 
  2629. | (after shifting).
  2630.     btst    IMM (0),d0    | is delta < 1?
  2631.     beq    2f        | if so, do not do anything
  2632.     tstl    d1        | is delta == 1?
  2633.     bne    1f        | if so round to even
  2634.     movel    d0,d1        | 
  2635.     andl    IMM (2),d1    | bit 1 is the last significant bit
  2636.     addl    d1,d0        | 
  2637.     bra    2f        | 
  2638. 1:    movel    IMM (1),d1    | else add 1 
  2639.     addl    d1,d0        |
  2640. | Shift right once (because we used bit #FLT_MANT_DIG!).
  2641. 2:    lsrl    IMM (1),d0        
  2642. | Now check again bit #FLT_MANT_DIG (rounding could have produced a
  2643. | 'fraction overflow' ...).
  2644.     btst    IMM (FLT_MANT_DIG),d0    
  2645.     beq    1f
  2646.     lsrl    IMM (1),d0
  2647.     addw    IMM (1),d2
  2648. 1:
  2649. | If bit #FLT_MANT_DIG-1 is clear we have a denormalized number, so we 
  2650. | have to put the exponent to zero and return a denormalized number.
  2651.     btst    IMM (FLT_MANT_DIG-1),d0
  2652.     beq    1f
  2653.     jmp    a0@
  2654. 1:    movel    IMM (0),d2
  2655.     jmp    a0@
  2656.  
  2657. Lround$to$zero:
  2658. Lround$to$plus:
  2659. Lround$to$minus:
  2660.     jmp    a0@
  2661. #endif /* L_float */
  2662.  
  2663. | gcc expects the routines __eqdf2, __nedf2, __gtdf2, __gedf2,
  2664. | __ledf2, __ltdf2 to all return the same value as a direct call to
  2665. | __cmpdf2 would.  In this implementation, each of these routines
  2666. | simply calls __cmpdf2.  It would be more efficient to give the
  2667. | __cmpdf2 routine several names, but separating them out will make it
  2668. | easier to write efficient versions of these routines someday.
  2669.  
  2670. #ifdef  L_eqdf2
  2671. LL0:
  2672.     .text
  2673.     .proc
  2674. |#PROC# 04
  2675.     LF18    =    4
  2676.     LS18    =    128
  2677.     LFF18    =    0
  2678.     LSS18    =    0
  2679.     LV18    =    0
  2680.     .text
  2681.     .globl    SYM (__eqdf2)
  2682. SYM (__eqdf2):
  2683. |#PROLOGUE# 0
  2684.     link    a6,IMM (0)
  2685. |#PROLOGUE# 1
  2686.     movl    a6@(20),sp@-
  2687.     movl    a6@(16),sp@-
  2688.     movl    a6@(12),sp@-
  2689.     movl    a6@(8),sp@-
  2690.     jbsr    SYM (__cmpdf2)
  2691. |#PROLOGUE# 2
  2692.     unlk    a6
  2693. |#PROLOGUE# 3
  2694.     rts
  2695. #endif /* L_eqdf2 */
  2696.  
  2697. #ifdef  L_nedf2
  2698. LL0:
  2699.     .text
  2700.     .proc
  2701. |#PROC# 04
  2702.     LF18    =    8
  2703.     LS18    =    132
  2704.     LFF18    =    0
  2705.     LSS18    =    0
  2706.     LV18    =    0
  2707.     .text
  2708.     .globl    SYM (__nedf2)
  2709. SYM (__nedf2):
  2710. |#PROLOGUE# 0
  2711.     link    a6,IMM (0)
  2712. |#PROLOGUE# 1
  2713.     movl    a6@(20),sp@-
  2714.     movl    a6@(16),sp@-
  2715.     movl    a6@(12),sp@-
  2716.     movl    a6@(8),sp@-
  2717.     jbsr    SYM (__cmpdf2)
  2718. |#PROLOGUE# 2
  2719.     unlk    a6
  2720. |#PROLOGUE# 3
  2721.     rts
  2722. #endif /* L_nedf2 */
  2723.  
  2724. #ifdef  L_gtdf2
  2725.     .text
  2726.     .proc
  2727. |#PROC# 04
  2728.     LF18    =    8
  2729.     LS18    =    132
  2730.     LFF18    =    0
  2731.     LSS18    =    0
  2732.     LV18    =    0
  2733.     .text
  2734.     .globl    SYM (__gtdf2)
  2735. SYM (__gtdf2):
  2736. |#PROLOGUE# 0
  2737.     link    a6,IMM (0)
  2738. |#PROLOGUE# 1
  2739.     movl    a6@(20),sp@-
  2740.     movl    a6@(16),sp@-
  2741.     movl    a6@(12),sp@-
  2742.     movl    a6@(8),sp@-
  2743.     jbsr    SYM (__cmpdf2)
  2744. |#PROLOGUE# 2
  2745.     unlk    a6
  2746. |#PROLOGUE# 3
  2747.     rts
  2748. #endif /* L_gtdf2 */
  2749.  
  2750. #ifdef  L_gedf2
  2751. LL0:
  2752.     .text
  2753.     .proc
  2754. |#PROC# 04
  2755.     LF18    =    8
  2756.     LS18    =    132
  2757.     LFF18    =    0
  2758.     LSS18    =    0
  2759.     LV18    =    0
  2760.     .text
  2761.     .globl    SYM (__gedf2)
  2762. SYM (__gedf2):
  2763. |#PROLOGUE# 0
  2764.     link    a6,IMM (0)
  2765. |#PROLOGUE# 1
  2766.     movl    a6@(20),sp@-
  2767.     movl    a6@(16),sp@-
  2768.     movl    a6@(12),sp@-
  2769.     movl    a6@(8),sp@-
  2770.     jbsr    SYM (__cmpdf2)
  2771. |#PROLOGUE# 2
  2772.     unlk    a6
  2773. |#PROLOGUE# 3
  2774.     rts
  2775. #endif /* L_gedf2 */
  2776.  
  2777. #ifdef  L_ltdf2
  2778. LL0:
  2779.     .text
  2780.     .proc
  2781. |#PROC# 04
  2782.     LF18    =    8
  2783.     LS18    =    132
  2784.     LFF18    =    0
  2785.     LSS18    =    0
  2786.     LV18    =    0
  2787.     .text
  2788.     .globl    SYM (__ltdf2)
  2789. SYM (__ltdf2):
  2790. |#PROLOGUE# 0
  2791.     link    a6,IMM (0)
  2792. |#PROLOGUE# 1
  2793.     movl    a6@(20),sp@-
  2794.     movl    a6@(16),sp@-
  2795.     movl    a6@(12),sp@-
  2796.     movl    a6@(8),sp@-
  2797.     jbsr    SYM (__cmpdf2)
  2798. |#PROLOGUE# 2
  2799.     unlk    a6
  2800. |#PROLOGUE# 3
  2801.     rts
  2802. #endif /* L_ltdf2 */
  2803.  
  2804. #ifdef  L_ledf2
  2805.     .text
  2806.     .proc
  2807. |#PROC# 04
  2808.     LF18    =    8
  2809.     LS18    =    132
  2810.     LFF18    =    0
  2811.     LSS18    =    0
  2812.     LV18    =    0
  2813.     .text
  2814.     .globl    SYM (__ledf2)
  2815. SYM (__ledf2):
  2816. |#PROLOGUE# 0
  2817.     link    a6,IMM (0)
  2818. |#PROLOGUE# 1
  2819.     movl    a6@(20),sp@-
  2820.     movl    a6@(16),sp@-
  2821.     movl    a6@(12),sp@-
  2822.     movl    a6@(8),sp@-
  2823.     jbsr    SYM (__cmpdf2)
  2824. |#PROLOGUE# 2
  2825.     unlk    a6
  2826. |#PROLOGUE# 3
  2827.     rts
  2828. #endif /* L_ledf2 */
  2829.  
  2830. | The comments above about __eqdf2, et. al., also apply to __eqsf2,
  2831. | et. al., except that the latter call __cmpsf2 rather than __cmpdf2.
  2832.  
  2833. #ifdef  L_eqsf2
  2834.     .text
  2835.     .proc
  2836. |#PROC# 04
  2837.     LF18    =    4
  2838.     LS18    =    128
  2839.     LFF18    =    0
  2840.     LSS18    =    0
  2841.     LV18    =    0
  2842.     .text
  2843.     .globl    SYM (__eqsf2)
  2844. SYM (__eqsf2):
  2845. |#PROLOGUE# 0
  2846.     link    a6,IMM (0)
  2847. |#PROLOGUE# 1
  2848.     movl    a6@(12),sp@-
  2849.     movl    a6@(8),sp@-
  2850.     jbsr    SYM (__cmpsf2)
  2851. |#PROLOGUE# 2
  2852.     unlk    a6
  2853. |#PROLOGUE# 3
  2854.     rts
  2855. #endif /* L_eqsf2 */
  2856.  
  2857. #ifdef  L_nesf2
  2858.     .text
  2859.     .proc
  2860. |#PROC# 04
  2861.     LF18    =    8
  2862.     LS18    =    132
  2863.     LFF18    =    0
  2864.     LSS18    =    0
  2865.     LV18    =    0
  2866.     .text
  2867.     .globl    SYM (__nesf2)
  2868. SYM (__nesf2):
  2869. |#PROLOGUE# 0
  2870.     link    a6,IMM (0)
  2871. |#PROLOGUE# 1
  2872.     movl    a6@(12),sp@-
  2873.     movl    a6@(8),sp@-
  2874.     jbsr    SYM (__cmpsf2)
  2875. |#PROLOGUE# 2
  2876.     unlk    a6
  2877. |#PROLOGUE# 3
  2878.     rts
  2879. #endif /* L_nesf2 */
  2880.  
  2881. #ifdef  L_gtsf2
  2882.     .text
  2883.     .proc
  2884. |#PROC# 04
  2885.     LF18    =    8
  2886.     LS18    =    132
  2887.     LFF18    =    0
  2888.     LSS18    =    0
  2889.     LV18    =    0
  2890.     .text
  2891.     .globl    SYM (__gtsf2)
  2892. SYM (__gtsf2):
  2893. |#PROLOGUE# 0
  2894.     link    a6,IMM (0)
  2895. |#PROLOGUE# 1
  2896.     movl    a6@(12),sp@-
  2897.     movl    a6@(8),sp@-
  2898.     jbsr    SYM (__cmpsf2)
  2899. |#PROLOGUE# 2
  2900.     unlk    a6
  2901. |#PROLOGUE# 3
  2902.     rts
  2903. #endif /* L_gtsf2 */
  2904.  
  2905. #ifdef  L_gesf2
  2906.     .text
  2907.     .proc
  2908. |#PROC# 04
  2909.     LF18    =    8
  2910.     LS18    =    132
  2911.     LFF18    =    0
  2912.     LSS18    =    0
  2913.     LV18    =    0
  2914.     .text
  2915.     .globl    SYM (__gesf2)
  2916. SYM (__gesf2):
  2917. |#PROLOGUE# 0
  2918.     link    a6,IMM (0)
  2919. |#PROLOGUE# 1
  2920.     movl    a6@(12),sp@-
  2921.     movl    a6@(8),sp@-
  2922.     jbsr    SYM (__cmpsf2)
  2923. |#PROLOGUE# 2
  2924.     unlk    a6
  2925. |#PROLOGUE# 3
  2926.     rts
  2927. #endif /* L_gesf2 */
  2928.  
  2929. #ifdef  L_ltsf2
  2930.     .text
  2931.     .proc
  2932. |#PROC# 04
  2933.     LF18    =    8
  2934.     LS18    =    132
  2935.     LFF18    =    0
  2936.     LSS18    =    0
  2937.     LV18    =    0
  2938.     .text
  2939.     .globl    SYM (__ltsf2)
  2940. SYM (__ltsf2):
  2941. |#PROLOGUE# 0
  2942.     link    a6,IMM (0)
  2943. |#PROLOGUE# 1
  2944.     movl    a6@(12),sp@-
  2945.     movl    a6@(8),sp@-
  2946.     jbsr    SYM (__cmpsf2)
  2947. |#PROLOGUE# 2
  2948.     unlk    a6
  2949. |#PROLOGUE# 3
  2950.     rts
  2951. #endif /* L_ltsf2 */
  2952.  
  2953. #ifdef  L_lesf2
  2954.     .text
  2955.     .proc
  2956. |#PROC# 04
  2957.     LF18    =    8
  2958.     LS18    =    132
  2959.     LFF18    =    0
  2960.     LSS18    =    0
  2961.     LV18    =    0
  2962.     .text
  2963.     .globl    SYM (__lesf2)
  2964. SYM (__lesf2):
  2965. |#PROLOGUE# 0
  2966.     link    a6,IMM (0)
  2967. |#PROLOGUE# 1
  2968.     movl    a6@(12),sp@-
  2969.     movl    a6@(8),sp@-
  2970.     jbsr    SYM (__cmpsf2)
  2971. |#PROLOGUE# 2
  2972.     unlk    a6
  2973. |#PROLOGUE# 3
  2974.     rts
  2975. #endif /* L_lesf2 */
  2976.  
  2977.